X-643-70-330 


NASA TM X-; 


INTERPLANETARY TRAJECTORY 
ENCKE METHOD 
FORTRAN PROGRAM MANUAL 
FOR THE I.B.M. SYSTEM/360 


SEPTEMBER 1970 



GSVC 


GODDARD SPACE FLIGHT CENTER 

GREENBELT, MARYLAND 

N70~ 


cs 

o 

<0 

s 

oc 

O 


(ACCESSJONNUMBER) 


MK2 a 

ItJMBFPl 


(PAGES) 


^ (NASA 


SA CR OR TMX OR AO NUMBER) 


(THRU) 

/ 


(CODE) 

3^ 


(CATEGORY) 



X-643-70-330 


INTERPLANETARY TRAJECTORY ENCKE METHOD 
FORTRAN PROGRAM MANUAL 
FOR THE I.B.M. SYSTEM/360 


September 1970 


GODDARD SPACE FLIGHT CENTER 
Greenbelt, Maryland 



ITEM PROGRAM MANUAL 
FORTRAN IV VERSION 


by 


Fred H. Whitlock 
Goddard Si)ace Flight Center 
Greenbelt, Maryland 

and 

Henry Wolf 
Leon Lefton 
Norman Levine 

Analytical Mechanics Associates, Inc. 

Jericho, New York 

A program of this nature must, of necessity, take a period of several 
years for its development; it is thus impossible to mention the names of 
all those who have contributed to its growth. It was originally conceived 
by S. Pines and H. Wolf at Republic Aviation Corporation under contract to 
NASA (NASA- 109) beginning in 1959* version is issued under contract 

to the Theoretical Mechanics Branch, Laboratory for Space Physics of the 
Goddard Space Flight Center. Major contributors have been C. Bergren, 

C. HLpkins, M. Wachman, L. Lefton, F- Shaffer, F. Whitlock and N. Levine. 

Numerous additions and improve iicj; ;s have been made to the current 
version including reprogramming fc;r •• le I.B.M. Operating System 36O 
Coii5)uters, and developnent is a cvmatlnuing effort. This edition covers 
the conversion of the original machins language program to Fortran IV. 

The program 2ias been and is available for general use to interested 
organizations.** The authors express their appreciat on to Mrs. Beatrice 
Boccucci: for her assistance in the : inal preparation of this manual. 


**Address: Head, Theoretical Mechanics Branch 

Laboratory for Space Physics, Code 6k3 
Goddard Space Plight Center 
Greenbelt, Maryland 20771 



TABI£ of comtemts 

Section 


I. Introduction I-l 

II. Notation II-l 

III. General Procedure for Using Program III-l 

IV. Initial Conditions IV-1 

A. Cartesian Coordinates IV-1 

B. Geodetic Polar Coordinates IV-2 

C. Geocentric Polar Coordinates IV-2 

D. Osculating Element Input IV-3 

E. Ccxmnents IV-3 

V. Terminating Conditions V-1 

VI . Permissible Perturbations VI-1 

VII. Badar Information Program VII-1 

VIII. Subroutine MODIF VIII-1 

A. Radiation Pressure VIII-1 

B. Aerodynamic Drag VIII-2 

C. Atmospheric Tables VIII-3 

D. Printout VIII-4 

E. Ephemeris Time VIII-5 

F. Inclusion or Exclusion of Perturbations VIII-6 

IX . Input IX-1 

X. Output X-1 

A. Program Outputs X-1 

B. Optional Outputs X-5 

1. Coordinates of Vehicle X-5 

2. Shadow Print . X-7 

3 . Radar Output X-8 

4. Reentry X -9 

5* Trajectory Search X -9 

6. Impact Parameters X-U 

7 . Apogee, Perigee, Nodal Crossing Print. . . . X-12 



TABLE OF CONTENTS (ContJ 

Section Psige 


XI. Internal Procedures XI-1 

A. Units XI-1 

B. Ephemeris Tape XI-1 

C. Ephemeris in Core XI-5 

D. Perturbation Program XI-4 

XII. Beferences XII-1 

XIII. Routines Available Upon Request XIII-1 

XIV. Methods of Integration XIV-1 

Appendices 

A. Introduction A-1 

B. Equations of Motion B-1 

C. Integration Procedure C-1 

D. Solution of the Kepler Two- Body Problem B-1 

E. Computation of Perturbation Terms E-1 

F. Conclusions F-1 

G. Oblateness Terms G-1 

H. Transformation Equations from Geodetic Polar Coor- 

dinates to Cartesian Coordinates H-1 

I. Transformation Equations for Radar Simulation . . . 

J. Variable Order Interpolation Scheme 

K. Drag Computation K-1 

L. Ccmputation of Subsatellite Point L-1 

M. Change of Independent Variable - Beta Mode M-1 

N. Shadow Logic N-1 

O. Solar Radiation Pressure 0-1 

P. Ecliptic Coordinates P”1 

Q. Moon Rotating and Fixed Coordinate System 

R. Trajectory Search . R-1 

S. Equations for Flight Path Azimuth and Flight Path 

Angle S-1 

T. Osculating Elements T-1 

U. Impact Parameters U-1 

V. Moon's Orbital Plane Input and Output V-1 

W. Equations for Translunar Plane Input W-1 



LIST OF FIGURES 

Page 

Figure 1. Geometry of the Elliptic Two- Body Orhit D-3 

Figure 2. Relation between Declination, Geocentric and Geodetic 

Latitudes H-2 

Figure 3 . Diagram for Shadow Logic N-3 



I . IMTROrUCTION 

0!his report describes a general purpose Interplanetary Trajectory 
Encke Method (ITEM) Program, programmed in the FORTRAN IV language. 

The method employed is designed to give the maximum available accuracy 
without incurring prohibitive penalties in machine time . On the basis of 
research described in Reference I 4 ., the Encke method was selected as best 
satisfying these requirements. However, the classical Encke method was 
modified to eliminate sane of its objectionable features. This modified 
Encke method is described in Appendix A. 

The perturbations included in this program are the gravitational 
attractions of the Earth, Moon, Sun, Mercury and the outer planets . The 
outer planets are considered as point masses. Additionally, the effects 
of the zonal and tesseral harmonics of the Earth, as well as aerodynamic 
dr€ig, small corrective thrusts, and radiation pressure including the 
shadow effect of the Earth, are considered. The input may be prepared 
in any one of several common systems and a great variety of output 
options are available . 

Additional options that are currently under development are explained 
in Section XIII. 
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II. NOTATION 

Upper case - vectors; Hats - unit vectors; Lower case - magnitudes 


Description 

Symbol 

Units 

Cartesian coordinatesof vehicle with respect to 
reference body 

X y z 

km 

Velocity components in Cartesian Coordinates 

• • • 

X y z 

km/sec 

Time 

t 

hrs . 

Longitude measured frcan Greenwich, + East 
(used in Section Iv and Appendix h) 

e 

degrees 

Longitude of vernal equinox 

00 

degrees 

Speed 

V 

km/sec 

Geodetic altitude* 

h 

km 

Geodetic latitude 

9 

degrees 

Geodetic flight path angle 

Y 

degrees 

Geodetic flight path azimuth 

A 

degrees 

Acceleration parameter (defined in Appendix E) 

u 


Right ascension 

RA 

degrees 

Astronomical units 

AU 


Earth radii 

ER 


Earth mass 

“e 


Semi -major axis 

a 

ER 


*Note: The following 5 symbols with primes denote the corresponding 

geocentric quantities. Geocentric in this report refers to a 
spherical earth, i.e., e2 = o. In this case cp’ 3 6 = 
declination. 
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Description Symbol 

Vehicle position vector R 

Distance to vehicle r 

Perturbation displacement vector AR 

Perturbation displacement vector components 

Perturbation acceleration F 

• . 

Coordinate functions and their time derivatives 
Mass parameter ji 

Earth's eccentricity as used in Appendices H, I, L, S e 

Mean motion n 

A A 

Unit vectors for classical two-body orbit solution P,Q 

Eccentric anomaly as used in Appendix T E 

Elevation angle as used in Appendix I E 

Ho * Ro <io 

Inclination of orbital plane i 

Right ascension of the ascending node n 

Vector from the body to vehicle R^ 

Greenwich Hour Angle G.H.A. 

Argument of perigee ou 

Parameters which account for polar oblateness of the 

earth, defined in Appendix H c, s 

Right ascension of the station meridian RA 

s 

Range measured from observation station p 

Direction cosines measured in a topocentric coordinate 

system X , p, 

Declination 6 
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SUBSCRIPTS 


Vehicle 

ith perturbing body 

Quantity obtained from Keplerian solution of 
two-body problem 

Reference body as used in Appendix B 
Station 



Value at rectification time 

Corresponds to x, y, z components respectively 
Value at perigee 


Symbol 

V 

i 


k 


c 

s 



o 


n - 1, 2, 5 
P 
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III. GENERAL PROCEDURE FOR USING PROGRAMS 

Initial conditions, terminal conditions and print frequency, as 
well as other parameters controlli::g the flow of the program, are read 
as input. The computation of the trajectory then proceeds until one of 
the terminal conditions ( e .g . maximum time) has been reached or an 
error is encountered. At this time the program prints the reason for 
its termination and proceeds to the next case. When an end of file is 
enccintered on the input tape, control is transferred to the monitor. 
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IV. INITIAL CONDITI' ' S 

The Initial ccnditions necessary for the specification of a 
trajectory are: 

1. Initial position of the vehicle relative to the reference body. 

2. Initial velocity of the vehicle relative to the reference body. 

3 . Initial time of launch referenced to a base tinge. 

For specification of the intial conditions^ the reference systems and 
units shown below may be used. 


Cartesian Coordinates 


The coordinate system is defined as follows: 

1. The origin is at the center of the reference body. 

2. The x-axis is in the direction of the meeui equinox of 
December 31-0 of the year of launch. 

3 . The xy plane is the mean equatorial plane of the Eeurth. 

Initial position is given by the x, y, z coordinates of the 
vehicle. Initial velocity is given by the x, y, z cchb- 
ponents of the vehicle . Initial time of launch from base 

time^^^ (t) is also given. If the program is used in its 

(2) 

standard form, the units to be used for the above are: 

X, y, z - kilcMneters 

X, y, z - kilometers/second 

t - month, day, hours, minutes, and seconds from 
base time 

The year of launch must also be given. 

U) 

( 2 ) 


Hie base time is 0.0*^ UT Deceniber 3 I of the year previous to the 
year of launch. 

Scale factors are used to convert.the input units to the units used 
internally (ER or AU a” other set of units may be used 

by changing the? sc'. ' . s with the appropriate ID card as 

described in Seoi c 
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B Geodetic Polar Coordipates 


Initial position of the vehicle is given by: 

1. Geodetic latitude (cp) 

2. Longitude' ^( 0 ), measured from Greenwich 
5* Geodetic altitude (h) 

( 3) 

4 . Longitude of vernal equinox'-^' at initial time (0o) • This 
quantity may be computed by the program or may be loaded 

Initial velocity of the vehicle is given by: 

1. Speed (v) with respect to the center of the Earth. 

2. Flight path azimuth (A) measured clockwise from north in 

a plane normal to the geodetic altitude. 

5« Flight path angle (y) measured from a plane normal to the 
geodetic altitude . 

Initial time of launch from base time^'^^t) must also be given. 
The following units must be used with the above quantities: 

1. cp# and 0 q - degrees; h - kilometers 

2. A and y - degrees; v - kilcaneters/second 
3 * t - month, day, hours, minutes and seconds 

C Geocentric Polar Ccx)rdinates 


Ordinarily an input given in polar coordinates will be interpreted 
as described in the preceding Paragraph B. However, if NOPT(l) =* 
the program w-t:i interpret latitude as declination, height as distance 
above a spherical Earth of equatorial radius, and flight path angle and 
azimuth with reference to a plane normal to the radius vector. 

( 3 ) If the right ascension (RA) at initial time is known, it may be used 

in place of longitude ( 0 ) , The longitude of the vernal equinox (0o) 
is then se^ to zero. 

(4) See INPUT Section. 
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D* Osculating Element Input 

The osculating elements to be input are: 

Argument of perigee 
Longitude of ascending node 
Inclination 

Semi-major axis (in Earth radii) 

Eccentricity 

Tims of perigee, mean anomaly, eccentricity anomaly, 
or true anomaly (only one) 

The program converts the above to Cartesian coordinates and then 
continues normally. (See Section IX, ID = 10.) 

E. Comments 

1. The program computes in Cartesian coordinates. Units used 
internally in the con^jutation are: 

(a) -Position; Earth Radii (ER) or Astronomical Units (AU) . 

(b) Velocity: ER/HOUR or AU/hOUR 

(c) Time ; Hours 

(Earth Radii units are used in the Earth or Moon 
Reference. Astronomical Units are used in the 
Sun, Mercury, Venus or outer planet reference) . 

2. The user is restricted to Cartesian coordinates when launching 
from any body other than the Earth as directed in preceding 
versions. 

3 . Equations for converting the initial conditions from Polar 
to Cartesian coordinates are shown in Appendix H. 
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V. TERMniATIMG CONDITIOMS 

The se-f of conditions which will terminate a trajectory may be 
summarized as: * 

1. Maximum time of flight - hours. 

2. Maximum distance from euiy possible reference body considered 
in the solution. last value in R- vector of integration and 
/rint block "ARRAY (l,n,j) ." 

3 . Minimum distance from any possible reference body considered 
in the solution. First value in R-vector of integration and 
print block ARRAY (l,n,l)* 

Any of these conditions will terminate a trajectory. Loading a 
large number into any of the maxima and a zero into any of the minima 
will make the corresponding conditions inoperative. A proper choice of 
these numbers will permit complete computation of the desired trajectory, 
avoid extensive unnecessary computation and guard against faulty input. 


* n designates reference 

n » 1 Jfcfeirth 
n - 2 Moon 
n => 5 Suii 
n » 4 Venus 
n = 5 Mars 
n = 6 Jupiter 


body. 

n ■ 7 Saturn 

n = 8 Uranus 

n = 9 Neptune 

n = 10 Pluto 
n = 12 Mercury 
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VI. PERMISSIBLE PERTURBATIOMS 

The trajectory computation consists of two parts, the exact solu- 
tion to the two-body problem and integrated additions to this solution 
for the effect of perturbations. The successful control of round-off 
errors in the modified Encke method depends on preventing the accumulated 
round-off error in the integrated perturbation displacement from affecting 
the computed position. This is achieved by always keeping the perturba- 
tion displacement small and rectifying whenever the perturbation exceeds 
specified limits. The constants mentioned below are used in determining 
the allowable limits as ratios of the perturbation position and velocity 
to the two-body position and velocity, respectively. 

This ratio for the position vector is shown in the following sketch. 



Trajectory 

Perturbation Displacement 
Kepler Trajectory 
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The recommended values for these ratios are as follows: 

Position Ratio 

POSRCS S .0001 

Velocity Ratio 

.2 

VELRCS < .0001 

and these are incc^orated into the program. Modifications may be made 
by altering the data subroutine, or by reading them in under ID = 12, 
or by using Subroutine Modif . 
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VII. RADAE INFORMATION PROGRAM 

The program may be used to simulate radar data if desired. A 
maximum of 50 stations can be processed at one time. The following 
information is required for each station considered: 

1. Station Name - for identification purposes 

2 . Position of Radar Station 

a. Longitude ( 0 ) of the station frcm Greenwich - degrees, 

minutes and seconds* - positive eastward. 

b. Geodetic latitude (cp) of the station - degrees, minutes 

and seconds* - positive north. 

c . Altitude ( h) of station above sea level - feet . 

The simulated radar information consists of azimuth, elevation, 
topocentric right ascension and declination, slant range, and range rate. 
It is printed at every normal print time for which the elevation angle 
is positive. Refraction is not considered. 

This section is coded as a subroutine and may be called at any time. 


* Alternatively, these quEuatities may be given in degrees and decimals . 
Zero's must be loaded into the positions reserved for minutes and seconds. 

The fractional parts will not appear in the print ouL reproducing the 
station coordinates. They will, however, be included in the computation. 
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VIII. SUBROUTINE MODIF 

Madifi cations to program constants, which normally remain unchanged 
during the running of a number of cases, may be made by using the 
appropriate common in conjunction with a compilation of a subroutine 
called MDDIF. This subroutine may include data statements, FORTRAN G or 
H level statements, and/or read statements. The use of read statements 
is suggested to facilitate stacking of cases . Any modification included 
in this subroutine will be operative for all succeeding cases, unless it 
is revoked. 

Modifications required more frequently may be accomplished through 
the use of ID(l) 's, as described in the INPUT section (Section IX) . 


A. 


Radiation Pressure may be included by loading a coefficient into 

RACOE ** 


The number to be loaded is; 

KC A 
r 

m 

is the radiation pressure in dynes/cm^ at a distance of 1 AU from the 
Sun. 

C =* 4.6 X 10 ^ .^2^— (estimated value) 

^ cm 


** Real number 
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A area in cm^ 

m mass In grams 

K scaler 5600^( 25455- )^/65?8. 165 X 10^ = .11178 x 10® seconds 
to hours, ER to AU, cm to ER 

The radiation pressure will only be active if sunlight impinges 
on the vehicle . For correct results the radiation pressure should 
therefore, be run only in conjunction with the optional shadow computa- 
tion as described in Appendix 0. 

If, however, the expected trajectory may be safely assumed to be 
entirely out of the Earth’s shadow, shadow testing may be avoided with 
a consequent saving in machine time. In this case, the following 
modification cards should be included in MODIF; 

SHTH = l.DO*^ 

SEMI » 1 .DO** 

nopt(iy) = 0 


B. Aerodynamic Drag 

If inclusion of the aerodynamic drag is desired, the drag parameter 
l/2 Cjj a/m may be initialized by ID * 23, as described in (Section IX) 
or loaded into subroutine MODIF by means of the following card: 

COEFL » ** 


The units of a/m are the area in cm^ and the mass in grams . A 
layered atmosphere rotating with the Earth is assumed. The density is 
obtained by a linear interpolation of the density-altitude table . The 
above may be Incorporated into the Block DATA subroutine. 


** Real number(s) 
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C . Atmospheric Tables 


Atmospheric tables for the drag computation are stored in core. 
They correspond to COSPAR International Reference Atmosphere of I96I, 
contained in the Report of the Preparatory Group for an International 
Reference Atmosphere, accepted at the COSPAR meeting in Florence, April 
1961. The units for the air density are grams/cm^ and the height is 
given in ER from the center of the Earth. If it is desired to change 
this atmosphere, the following procedure should be followed: 

NTAR ■ * The number to be entered 

is N - the number of 
entries in the density table. 


RTBUI) =» ** I » 1, 2, - - -, N - the values 

of r for which densities are 
given, in ascending order (a 
maximum of 50) . 


RHOT(i) • ** I = 1, 2, , W - the values 

for the air density in grams/cm® 
in respective order corresponding 
to the preceding r table . 


If other units are used for the density table, the drag parameter des- 
cribed in Pajrt B of this section must be read in with like units and the 
constant ( -6578.165D5)** normally in DRSC has to be changed accordingly. 
The negative sign directs the drag force opposite to the velocity. This 
constant converts the drag from the units used for A, M and p to ER/hour^. 


* Integer 
** Real number 
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D. Printout 


Ihe program provides a special printout near the Earth, Moon, Sun, 
or T- Planet reference. This prinout occurs at every integration step 
and is useful for observing the behavior of these relevant quantities 
during ascent and re-entry. This feature is triggered by the following 
modification cards: 


SERE(l) = 

** 

I » 1, 2, . . 

For printout near Earth use 

index 1 

.*a ER units 

Me on 

2 

ER 

Sun 

3 

AU 

Venus 

4 

AU 

Mars 

5 

AU 

Jupiter 

6 

AU 

Saturn 

7 

AU 

Uranus 

8 

AU 

Neptune 

9 

AU 

Pluto 

10 

AU 

Mercury 

12 

AU 


12 


The numbers used above are the radial distances within which the special 
printout is to be effective. The uoits are earth radii for the Earth and 
Moon references and astroncmilcal units for the remaining planet references. 
A zey'O in any of the SEPE(l) *s suppresses this feature. 


** Real number(s) 
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E. Ephemerls Time 

The planetary coordinates are interpolated using ephemeris time. 

ET = UT + A T E 

An approximate value of A T E (35 seconds) is used. To change this 
quantity, the following card, giving A T E in hours, is Inserted: 

ETMUT = ** The value for Ate- hours. 

To restore original quantity: 

ETMUT = .009888888889^*^ A T E is 35 seconds. 


** Real number( s) 
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F. Inclusion or Exclusion of Perturbations 


Ordinarily included are the graviational attractions of the Moon, 

Sun, Mercury, and the outer planets. The gravitational field of 
the earth (2nd, 3^^ 4th zonal harmonics) are Included if 
M0PT(50) =1.* A maximum of 15 zonal coefficients may be used. To 
exclude any or all of these perturbations, the following modifications 
should be included in subroutine MDDIF: 



MOON 

MEI(2) = O.DO 


SUN 

MEI(3) =» O.DO 


VENUS 

MEI(4) » O.DO 


MASS 

MEI(5) » O.DO 


JUPITER 

MEI(6) « O.DO 


SATURN 

MEI(t) =* O.DO 


URANUS 

MEl(8) » O.DO 


NEPOUNE 

MEI(9) -O.DO 


PLUTO 

MEI(IO) » O.DO 


MEaRCURY 

MEI(12) = O.DO 


ZONAL and/or TESSERAL HARMONICS** 

Zonals 

EFJ 

..15) = O.DO 

Tesserals 

ECG(I,J) » O.DO 

■pOn/ T /> T\/\ 

^ I and J ■ 1 to 


ESG(I,J) = O.DO 


* Described in the INPUT Section (Section IX) 
** Described in Appendix G 
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IX. mpuT 

Input to the program is read in as follows: 

Each set of input is preceded by an ID card which contains an integer 

number terminating in column 10. This card may also contain Hollerith 

information start ir^ in column 11. 

ID » 1 Permits one card of Hollerith information - usually used for 
case identification. 

ID * 2 Permits one card containing a set of 72 fixed point I's or 
O’s. Each flag (l or o) corresponds to the same numbered 
subroutine. A zero is used for normal operation and a one 
is used to print diagnostic information in the proper sub- 
routine. A blank card after ID = 2 will be necessary if the 
system does not zero out core before load time and normal 
operation is desired. In the program, these flags are 
referred to eis NC( l) . 

ID » 3 Performs similarly to ID = 2. It allows a card of I's and 
O's to be read into NOPT(i) (I = 1 to 72) . Oliese flags 
permit the incorporation of various options into the 
program. Following are the currently available options: 
N0 PT(i) =* 1 indicates polar geocentric coordinates 

» 2 indicates geodetic coordinates when polar 
load is used. 

NOPT ( 2 - 15 ) are used for print options. 

* 0 indicates no print 
» 1 indicates print 

N0PT(2) is associated with XR print 

N0PT(5) is associated with XRDT print 

N0PT(4) is associated with XVE print 

N0PT(5) Is associated with XVM print 

N0FT(6) is associated with XME print 
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N0PT(7) is associated with XVS print 

N0PT(8) is associated with XVYN print 

N0PT(9) is associated with XVMR print 

NOPT(IO) is associated witn XVJP print 

NOPT(H) is associated with XI print 

N0PT(12) is associated with XIDT print 

N0PT(13) is associated with lEXI print 

N0PT(14) = 1 prints data statement parameters 

N0PT(15) = 1 deletes regular print after rectification 

N0PT(i 6) » 1 deletes print in rectification 

N0PT(17) = 1 activates shadow computations 

N0PT(20) a 2 activates predictor only integrator* 

N0FT(20) » 3 activates predictor-corrector integrator* 

N0PT(25) = 1 prints spin axis angles. This option may he 

used in ccnjunction with NC(l<-7) to allow the user 
to input a desired spin vector. 

N0FT(33) = 2 rotates positions and velocities from equatorial 
to ecliptic coordinates for printing 
N0PT( 40-44) = 0 indicates no print 

N0PT(40 and 4l)= 1 is associated with XVPl print 

N0PT(40 and 42) =1 is associated with XVP2 

N0PT(40 and 43) *1 is associated with XVP3 

N0PT(40 and 44) »1 is associated with XVP4 

N0PT(40 and 45) ■! is associated with XVP5 

N0PT(40 and 46) =1 is associated with XVP6 

N0PT(50) » 0 No zonal or Tesseral harmonics are used. 

N0PP(50) = 1 Includes the gravitational field of the Earth 
(2nd, 3rd, and 4th zonal harmonics) 

N0PT(50) = 2 Includes zonal and tesseral harmonics 


* Described in Section XIV Methods of Integration 
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N0PT(51) = 1 
N0PT(65) =* 1 
N0PT(68) = 1 

N0PT(68) = 2 

N0PT(69) = 1 

N0PT(70) = 1 


IELD=2 


IELI>=3 


NOPT(TO) =» 2 
N0PT(70) = 3 

N0FT(71) = 1 
N0PT(72) = 1 


activates restart feature (See ID»9) . 
activates trajectory search (lD=l8 ipust be included) . 
activates residual computations with radar stations. 
(ID=21 must be Included) . 

activates residual computations without radar stations. 

(ll>>21 must be included) . 
activates solar engine logic for small corrective 
thrusts. QSiis option can be made available 
upon request. 

activates element roatico. This option is used in 
conjunction with a variable called lELD. lETiT) 
should be set in subroutine MODIF as follows: 
Permits ecliptic elements for input, resulting in 
equatorial osculating elements outjHit. 

Permits equatorial elements for input, resulting in 
ecliptic element output. 

Permits osculating element input, Brouwer mean 
element output. 

This option cein be made available upon request, 
activates Beta Integrator* 

permits element load with the period (HRS) sub- 
stituted for the semi-major eucis A(ER) . 


* Described in Section XIV Methods of Integration. 



ID = 4 


Used to read in start time of flight in year, month, day, 
hours, minutes, and seconds: starting reference body; and 

target reference using the following format (515> F5*2, 215). 
Bie reference bodies are numbered as follows: 

1 = Earth 4 = Venus 7 = Saturn 10 = Pluto 

2 = Moon 5 = Mars 8 = Uranus 12 = Mercury 

5 = Sun 6 = Jupiter 9 = Neptune 

ID =» 5 Used for polar load and reads in 0, cp, h, v. A, y, and 0 q. 

The format used is (3D20.0) . The progreun expects all 
angles in degrees, altitude in kilometers measured from 
the surface of the Earth, and the velocity in kilometers 
per second. If 0o is read in as 1000. ODO, the program 
Vill compute the proper Gq. 

ID ■ 6 Used for Cartesian ?nput. x, y, z, x, y, and z are read 
in with format (3E20.0) . The program expects these coor- 
dinates to be equatorial in kilometers and kilometers per 
second, with the starting reference body as center. 

ID = 7 This option generates initial conditions for a trajectory 
which is designed to get a spacecraft to the target within 
a specified number of dsys, without thrust. The input is: 
the Julian date of start time, the flight time in days, 
option number, and the radius of the Earth's sphere of 
influence. The format used is (3D20.0) . 

Option number 1 starts in Sun Reference . 

Option number 2 and 3 are not activated, however, 
they can be made available upon request. (See 
Section Xm) . 



When Option » 1 is used, input for II>=7 “ust be followed by 
4 cards 

Card 1 =» 0, 0, 0 .51373647D-6* Format (l3,D12.6,Dl8.8) 
Card 2 ^ Blank 
Card 3 =* Blank 
Card k = Blank 

ID = 8 This ID permits one to read in a vector of special print 

times . The first card after the ID card contains the num- 
ber of such print times from 1 to 50, format (llO) . 

Qlie following card or cards contain the times, format 
(3De0.0) . If this ID is used in conjunction with N0PT(69) 

( solar engine option) , these times are used forstarting 
and stopping the solar engine . Odd mambers start the 
engine, the even numibers shut the engine off. 

ID = 9 This ID reads in the following: PRSP, RACOE, TIMED, STI, 

VI, CCHT, ENPLAN, and TR using format ( 5 DBO.O) . 

PRSP A non-zero value will suppress normal print times 

until the specified time in hours has been reached. 

RACOE A non- zero value will activate radiation pressure 
canputations . 

TIMED Maximum time of flight in hours. 

STI Not activated 

VI Not activated 

CCNT Triggers nodal crossing print. (See ID = l4 also). 

If CCNT and ICCNT are set to non-zero, ICCNT will 
take precedent . / * 

ENPDAN A non-zero activates the number of perturbing bodies 
to be used in the calculations. This must be an 
integer 1 or 12. If 1 is used the ephemeras data 
tape is neigher required nor used. 

Gravitational constant for the Sun 
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TR A non- zero value, in conjunction with N0PT(51)^ will 
activate the restart feature. TR should be set to 
the desired restart time in hours. 


ID = 10 This ID permits one to load the initial conditions as 

osculating elements of an ellipse . The following are read 
in: SOMEG, LOMEG, INC, A, ECC, ELOAD, ELTRIG with format 

(3D20.0) . 

SOMEG is the argument of perigee 

LOMEG is the longitude of the ascending node 

A is the semi-major axis 

ECC is the eccentricity 

ELOAD depends on ELTRIG 


If ELTRIG * 1 
ELTRIG * 2 
ELTRIG = 3 
ELTRIG » 4 


ELOAD = time of perigee 
ELOAD = mean anomaly 
ELOAD = eccentric ancxnaly 
ELOAD = true anomaly 


ID =» 11 Permits one to alter the integration and print intervals of 
the various reference bodies. The number of cards to be 
read is a function of ENPLAN. (See Sample Input Data in 
Section IX ) . The data expected are read in with format 
(3D20.0), and terminated by any integer greater than or 
equal to 13 - (format IlO) as follows; 

Card 1 contains eight distances from Barth in ER 

2 contains seven integration intervals in hrs. 

3 contains seven print intervals in hrs. 

4 contains eight distemces from the Moon in ER 

5 contains seven integration intervals in hrs . 

6 contains seven print intervals in hrs. 

7 contains eight distances from the Sun in AU 

8 contains seven integration intervals in hrs . 

Q contains seven print Intervals in hrs. 
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10 contains eight distances from Venus in AU 

11 contains seven integration intervals in hrs. 

12 contains seven print intervals in hrs. 

15 contains eight distances from Mars in AU 

lif contains seven integration intervals in hrs. 

15 contains seven print intervals in hrs . 

16 contains eight distances frcM Jupiter in AU 

17 contains seven integration intervals in hrs . 

18 contains seven print intervals in hrs. 

LAST contains an integer .GE. I5 (FORMAT IlO) . 

ID » 12 Permits one to make changes in the program's built-in data 
or to read in other-than-normal inputs. This can be done 
by using a subroutine called I^DIF which must contain the 
proper block ccMmion. 

ID = 13 Allows one to change input and output scale factors, using 
format (5D20.0) . The card following the ID card contains 
TSCL, HEKM, and XMEfCM. 

TSCL is the time scale factor and sits in the program 
as 5600. It is used to change seconds to hours 
and hours to seconds 

REKM sits in the program as 6378 *165^ the number of 
kilcxoeters in one ER. 

XMEKM sits as 14*9599 X 10^ and is the number of kilo- 
meters in one AU. 

ID = 14 sets triggers for apogee, perigee, and nodal crossing prints 
The following card reads in ICART, ICPNT, and ICCMT, using 
format (3II0) 

*fcll 

ICANT n prints every n^ apogee 

'til 

ICPNT = n prints every n perigee 

til 

ICCHT = n prints every n^ nodal crossing 
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ID » 15 Not activated 

ID ■ l6 Allows radar station data to be read in. The card following 
the ID card contains the number of stations to be read in 
with format ( IlO) . !Che next two cards contain the name eind 
coordinates of the first station with format (A46/7niO.O) . 
This last fonnat is repeated until all stations have been 
accounted for. 

ID = 17 is used for reading in solar engine information with cor- 
corrective thrusts. This option is not activated*. 

ID = 18 is used for the iterator. The card following the ID card 
is read in with the fonnat (IOI5) 

IPS =« Number of the dependent parameter(s) . 

UMAX = Maximum number of iterations. 

NSL «» Number of the independent parameter( s) . 

IPOFiXiPS) ■ Identification number(s) of the quantities to 
be achieved. 

IVAfi(NSL) » Input quantities to be varied in the units and 
sequence of the polar load, i.e., 6,cp,h,v,A,y 

The remining cards are read in with the format (7DIO.O) 
YEPS(lPS) = The tolerance to which convergence is desired. 

if the solution converges to within the special 
tolerance, the iteration will stop. 


* See Section XIII 
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xeps(nsl) = 
xvar(nsl) « 

YC0N(IPS) «• 


Epsilon values of the independent va iable . 
Epsilons to he used for generating secant 
partials . 

The desired values of the dependent variables . 


A maximum of any six dependent variables may be selected. The identi- 
fication number(s) of the quantities to be achieved are as follows: 


1. B • T 

2. B • R 

3 . Earth, 

4 . Earth, 

5. Earth, 

6. Earth, 


Miss parameters 
Miss parameters 

Lunar or’ T- planet inclination.* 

Lunar or T- planet ascending node .* 
Lunar or T-planet argument of perigee.* 
Lunar or T-planet pericynthion radius.* 


ID » 20 Starts the program. 


ID = 21 Used in conjunction with N0PT(68) for reading in nominal 
and variational triggers with (exit) radar simulation. 

If N0PT(68) = 1, the card following the ID card is read 
with format (4ll0) . 

NOMTRI 1 = Nominal trajectory, 0 = variational 
NOQAN » The number of radar stations 

lOSCTR 0 =« Radar residuals, 1 = osculating element residuals. 


Qlie user may find this option convenient for checking out 
numerical derivatives. If N0PT(68) ■ 2, the card following 
the ID card is read with format ( UO) . 

NOMTOI 1 Nominal trajectory, 0 = variational 


* See Appendix R. 



ID = 25 Used to initialize aerodynamic drag computation 

The card following the ID card is read in with format 
(5D20.0) . 

COEFL* 1/2 a/m. 

DRSCf*^ 

DNTAR*^ number of entries in the density table. 

ID =» 4 must precede ID »= 5* Except for this condition and 
ID * 20 which must be read in last, the ID's may be read 
in randomly. 


♦ See Section VIII- B 
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100 c 

200 C 
300 C 
10100 
10200 

10300 

10301 

10500 C* 
10600 
10700 
10800 
10900 
11000 


20100 
202C0 
20300 
20400 
20500 C* 
20600 
20700 
20800 
20900 
21000 


SAMPLE LIST OF SUBROUTINE MODIF 


SUBROUTINE MODIF 
IMPLICIT REAL *8 (A-H, O-Z) 

COMMON/HENRY/H 1 (480) , MEI (1 2) , H2(226) 

REAL *8 MEI 

A** *********w**-*A«r******-t***<lr*'*********#***#*'*'*************#****** 


MEI (4) = O.ODO 
MEI (6) * O.ODO 
MEI (9) = O.ODO 
RETURN 
END 


SUBROUTINE MODIF 
IMPLICIT REAL *8 (A-H, O-Z) 

COMMON/GLOBE/FJG(15),CG(15,15),SG(15,15), EF ECG(15,15), 

XFSG(15, 15), EJG{4), ES22, EC22, ES31 , EC31 , EabJ, EC33, NIZ, N2S , N3T,NOQ 
*** *************************************************************** 


EFJ(2) = 1082.30D-6 
EFJ(3) = -2.30D-6 
EFJ(4) = -1.80D-6 
RETURN 
END 
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SAMPLE INPUT DATA 


1//GO.DATA5 DD 

2 I 

3 

4 

5 

6 


TEST CASE FOR I 
NOPT(50) = 1 


B.M. 360 VERSION 
FOR OBLATENESS 


000000000000000000000000000000000000000000000000000000000000000000000000 

3 

211111111111 110000000000000000003000000000000000010000000000000000000000 


8 


6 

CARTESIAN LOAD. 


9 

6095.8926D0 

604.36567D0 

2406.3589D0 

10 

1 . 101985D0 

9.8589445D0 

-4.4549856D0 

11 


4 



12 

1965 

5 30 

12 70.0 1 4 


13 


14 



14 


0 

0 0 


15 


9 



16 

0.000 


O.ODO 

10. ODO 

17 

O.ODO 


0. ODO 

O.ODO 

18 

1 .ODO 


O.ODO 


19 


1 1 SET VEaORS FOR DESIREH PLANETS 


20 


1 

SET EARTH VECT 


21 

1 .ODO 


4. ODO 

8. ODO 

22 

125. ODO 




23 





24 

0.0625D0 

0.0625D0 

0.5D0 

25 





26 





27 

1 .DO 


2. DO 

6. DO 

28 





29 





30 


13 

SET K TO TERMINATE READ. 


31 


16 

RADAR STATIONS 


32 


2 



33 

ROSMAN 




34 

277 

.0 7.0 

4.2 35. 

12. .7 

35 

CANBERA 




36 

148. 

57. 

20.9 -35. 

37. 52.7 

37 


20 STARTS THE PROGRAM. 


38/ 

* 





2874.56 

3766.588 


/ 
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1 //GO.DAT.\5 DD • 

2 1 ITERATOR DATA FOR MARS TRAJ. 1973 

3 NOPT(65) = 1 FOR SEARCH. 

4 2 

5 00000000000000000000000000000000000000000000000000000000000000000000000 

6 3 POLAR GEOC 

7 1 1 110010100001110000000000000000000000000000000001000000000000001000000 

8 4 START TIME 



0.0 

1 . 
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X. OUTPUT 


A. Program Outputs 

The following information is printed as the output of the program. 

1. Title 

2. Case number and any identifying titles. 

5 . Launch time - year, month, day, hour, min, sec. 

4. Input - in the same units as they were entered into the program. 

5 . List of parameters used in run. 

6. At each rectification the following data are printed; 

(b) RECTIFICATION PRINT (a) REFERENCE 
PERT OVER UNPERT » ( c) TIME « (d) DELTA T = (c) 

( a) Reference body 

(b) and (c) indicate the reason for rectification 

(c) If (c) * 0, rectification may be due either to switch of 
reference body or to change of integration interval. 

If (c) ^ 0, then the position, velocity, perturbations or 
the incremental eccentric ancmialy have exceeded the permis- 
sible limits and (b) indicates which has been exceeded 
( see Section Vl) . These indications are given as ; 

PO - Position 
VL - Velocity 

TH - Incremental eccentric anomaly 

(d) Time of rectification 

(e) Integration interval 
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TIME IN YEAR, MONTH, Bkx, 

HOUR. MIN, SEC, 

(a) 

T = HOURS FRCM EPOCH 

(b) 

XR YR ZR 

RR 

(c) 

XRDT YRDT 

ZRDT RRDT 

(d) 

RIGHT ASCENSION ( DBG) • _ 

DEXHi * 

(e) 

SUBSAT POINT LONG 

m 

(f) 

LAT 

= 

(g) 

HT 

= 

(h) 

GHA 

3 

(i) 

GEOCENTRIC AZIMUTH 

= 

(j) 

ELEVATION 

= 

(k) 

GEOIffiTIC AZIMUTH 


U) 

ELEVATION 

= 

(m) 


(a) Year, Month, day, hour, second frcm time of launch 

(b) Print time in hours from time of launch 

( c) Position coordinates and magnitude of radius vector with respect 
to the reference body - kilaneters/second 

(d) Velocity components and magnitude of velocity vector with respect 
to the reference body - kilometers /second. 

(e) Bight ascention and declination in Earth reference system - degrees 

(f) Longitude or sub- satellite point - degrees 

(g) Latitude (geodetic) - degrees 

(h) Geodetic hei^t above the Earth's surface - kilometers 

(i) Greenwich hour angle - degrees 

( j) Geocentric flight path azimuth - degrees 

(k) Geocentric flight path angle - degrees 
(-t) Geodetic flight path azimuth - degrees 
(m) Geodetic flight path angle - degrees 
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MOON SUBSAT POINT 


LONG 
LAT 
AZIM 

ELEIV =» 

OSCULATING ELEMENTS AT TIME T » 

TRUE ANOMALY = 

SEM MAJ AXIS = 

ECGENT » 

PERICENT * 

APOCENT = 

INCLINATION « 


(a) 

(b) 

(c) 

(d) 


(e) 

(f) 

(g) 

(h) 

(i) 

(J) 


(a) Moon longitude - angle between the projection of the vector 
from the Moon to the vehicle onto the Moon's orbital plane 
and the Moon-Earth vector (licx)n reference only) - degrees 

(b) Moon latitude - angle between the radius vector connecting 
the Moon and the vehicle and its projection onto the orbital 
plane of the Moon about the Earth (Moon reference only) - 
degrees 

( c) Selenocentric flight path azimuth - degrees 

(d) Selenocentric flight path angle - degrees 

(e) True anomaly - degrees 

(f) Semi-major axis of trajectory - ER _ ^ hyperbola 

(g) Eccentricity of trajectory** 

(h) Closest distance to the reference body (not necessarily the 
Earth)** - kilometers 

(i) Farthest distance frcm the reference body (not necessarily the 
Earth) **( meaningful only for elliptic orbits) - kilometers 

( j) Inclination of the orbital plane defined as the angle between 
the positive polar axis and the angular mcmentum vector** - 
degrees 


** These are the osculating values and hence only constitute an estimate 
of the quantities described. 



ARG PERIC 

= 

(a) 

PERIOD 


(b) 

MEAN MDT 

= 

(c) 

R A ASC NODE 

=j 

(d) 

M ANOMALY 

= 

(e) 

E ANOMALY 

9 

(f) 

T PERIC 

9 

(g) 

UNIT lERICENTER POSITION VECTOR » 

(h) 

UNIT ANGULAR MOMENTUM VECTOR = 

(i) 


(a) Argument of pericenter - angle measured from the ascending 
node to the pericenter vector**- degrees . 

Set to zero for circular orbits and poorly determined for near- 
circular orbits. 

( b) Period**r hours 

(c) Mean motion** - radians/hour 

(d) Right ascension of the ascending node measured from the vernal 
equinox eastward along the equator**- degrees 

(e) Mean anomaly** - degrees 

(f) Eccentric anomaly** - degrees 

(g) Time of nearest pericentei** - hours 

(h) Components of the unit vector directed from refexence toward 
pericenter** 

(i) Components of the unit angular momentum vector 


** Osculating values 


X-4 



B. Optional Outputs 


1. XVE 

9 

YVE 

a 

ZVE 

EVE 


(a) 

XVM 


YVM 

M 

ZVM 

RVM 

a 

(b) 

XME 


YME 

M 

ZME s 

RME 


(c) 

XVS 


YVS 


ZVS 

_ RVS 


(d) 

XWN 

3 

YWN 


ZWN = 

RWN 

s 

(e) 

XVMR 

m 

YVMR 


ZVMR ■ 

RVMR 

m 

(f) 

XVJP 


YVJP 


ZVJP • 

R7JP 

a 

(g) 

XVPl 


YVPl 

m 

ZVPl * 

RVPl 

a 

(h) 

XVP2 

SB 

YVP2 


ZVRJ 

RVP2 


(i) 

XVP3 

m 

YVP3 

sz 

ZVP3 = 

R7P3 

a 

(J) 

XVP4 


YVP4 


ZVP4 ’ 

RVP4 

a 

(k) 

XVP5 


YVP5 


ZVP5 = 

RVP5 


U) 

XVP6 

_ 

xvp6 

s 

Z7P6 • 

RVP6 


(m) 

XI 


ETA 


ZETA = 

PERT 


(n) 

XIDT 


EPADT 


ZETADT = 

VPERT 

M 

(o) 

D2XI 

s 

_ reETA 

m 

D2ZETA = 

APERT 

S 

_ (p) 


The above optional output appears between XRDT and RIGHT ASCENSION 
in the standard output. For instructions on how to obtain, see 
Section IX, ID » 5* 

(a) Coordinates of vehicle with respect to the Earth - kilometers 

(b) Coordinates of vehicle with respect to the Moon - kilometers 

(c) Coordinates of the Moon with respect to the Earth - kilometers 

(d) Coordinates of vehicle with respect to the Sun - kilometers 

(e) Coordinates of vehicle with respect to Venus - kilometers 

(f) Coordinates of vehicle with respect to Mars - kilometers 

(g) Coordinates of vehicle with respect to Jupiter - kilometers 
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(h) Coordinates of vehicle with respect to Saturn - kilometers 

( i) Coordinates of vehicle with respect to Uranus - kilometers 

( j) Coordinates of vehicle with respect to Neptune - kilometers 

(k) Coordinates of vehicle with respect to Pluto - kilometers 

(-t) Coordinates of vehicle with respect to E-M Barycenter - kilometers 

(m) Coordinates of vehicle with respect to Mercury - kilometers 

(n) Perturhation vector and magnitude of the perturbations with 
respect to the reference body - kilcxneters 

( o) Perturbation velocity vector and magnitude - kilometers/second 

(p) Perturbation acceleration vector and magnitude - kilcxneters/second^ 
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2. Shadow Print* 


PASSAGE FROM 


SHADOW 

PENUMBRA 

SUN 

PENUMBRA 


TO 


PENUMBRA 

SHADOW 

PENUMBRA 

SUN 


SHADOW 

AT (a) TIME IN ^^UMBRA 

PENUMBRA 


(b) ACCUMULATED TIME (c) 


The above optional output appears before TIME IN YEAR, MONTH, DAY, HOUR, 
MIN, SEC in the standard output. It is controlled through the INPUT 
subroutine [N0PT(17)] (see Section IX, ID » 3) . 


(a) Time at which vehicle traverses denoted shadow boundary - hours 

(b) Total time the vehicle spends in denoted shadow region during 
current traverse - hours 

(c) Total accumulated time spend in denoted shadow region since 
launch - hours 


* Apper iix N 
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3 • Radar Outpuf ^ 


STATION 

(a) 

AZIMUTH 

(b) 

ELEVATION 

(b) 

TOPOC. R A 

(c) 

TOPOC. DECL. 

(c) 

SLT RNG 

(d) 

RANGE 

(e) 


This output appears at the tail end of a normal printout . An ID 
card in the INPUT subroutine will control this segment of the 
program (see Section IX, ID =» l6) . 

(a) Station name (identification) for each station 

(b) Azimuth and elevation with respect to each station - degrees 

(c) Topocentric right ascension and declination with respect to 
each station - degrees 

(d) The slant range to each station - kilometers 

(e) Rate of change of slant range for each station - kilometers/ 
second 

If the elevation is negative (the vehicle is below the horizon), 
this print is suppressed for the station in question. 


* Appendix I . 
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4. Beentry Output 


REENTRY PRIWr TIME INERTIAL SPEED 

( kilometers / second) 

Right ascension, declination. Earth subsatellite points and flight 
path azimuth and angle as given above. 

The above optional output appears between GEOD ELE7 and MOON 
SUBSAT POINT in the standard output. 


5 . Trajectory Search Output 

The output consists of the normal ITEM output for a nominal trajectory 
and the same trajectory output for each variation requested for each 
iteration. The output format used only for the trajectory search 
follows : 

VARIATION IN INITIAL 

CONDITIONS (a) (b) ( c) (d) (e) (f) (g) 

(a) Change in latitude - degrees 

(b) Change in longitude - degrees 

( c) Change in altitude - kilometers 

(d) Change in velocity - kilometers /sec and 

(e) Change in azimuth - degrees 

(f) Cnange in flight path angle - degrees 

(g) Change in initial time - hours 
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5* Trajectory Search Output (cont.) 


QUANTITY CODE (a) 

DESIEED VALUES OF ABOVE 

QUANTITIES (b) 

REQUIBED ACCURACY (c) 

MATRIX OF PARTIAL DERIVATIVES (d) 

RESm^iUS AND HCANGES IN 

INITIAL CONDITIONS (e) (f) 



(a) Code indicating quantities to be searched for. 

(b) Desired values of above quantities - degrees, kilometers, 
seconds. 

(c) Tolerances allowed on above values - degrees, kilometers, 
seconds . 

(d) Jfetrix with the dependent variables arranged by row. The 
Independent by column. 

(ej Residuals (desired-nominal) of quantifies designated by the 
quantity code . 

(f) Change required in initial conditions. 

(g) Normalized changes in initial quantities in order of the 
variations . 


The option associated with trajectory search routines is initiated 
by an ID card in the INPUT subroutine ( see Section IX, ID = l8) 
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6 . Impect Pareuaeter Output * 


SHAT 

RHAT 

THAT 

E 

B*T, B’R 


(a) 

(b) 

(c) 

(e) (f) 


(a) Itoit vector in the direction of the incoming asymptote. 

(h) Oait vector normal to TEAT and the asymptote in a right-hand sense. 

(c) Unit vector parallel to the ecliptic (or Moon orbital plane) and 

perpendicular to the incoming asymptote. 

(d) Ve>, or from the Dody to the vehicle as it crosses the Impact plane. 

(e) The dot product of at the crossing and THAT. 

(f) The dot product of R at the crossing and RHAT. 


♦ . 'irdix U 
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7 * Apogee, Perigee, Nodal Crossing Print 

Apogee and perigee print times are ccaaputed more accurately 
than formerly, and the nodal crossings are found by iteration. 

The trigger for apogee is ICAHT; for perigee ICIWP; for nodal 
crossing ICCRT. The apogee and perigee print times are found by 
fitting a parabola through three neighboring points and deter- 
mining the miniimnn or maximum respectively. The crossing time is 
found iteratively by 


t ■ ® t 

CN CN-1 




t- is the time for which Z changed in sign. 
00 
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XI INTERNAL PROCEDURES 
A. Iftiits 

15ie units used internally are Earth Radii eind Barth Radii /hour 
in the Earth and Moon references, and Astronomical Iftiits and Astro- 
nomical Units/hour in the Sun, Mercury, Venus thru Pluto Reference 
systems . 

3. Ephemeris Tape 

The relative positions of the solar system bodies are obtained 
from a tape generated by the Jet Propulsion laboratory. A separate 
program prepares a binary tape referred to the mean equinox of MID- 
Flli^, containing l6 days perrecord, in a form compatible with the 
main program. The Ephemeris Subroutine searches the tape and 
reads in the proper file and record, keeping 52 days of tables in 
core storage at a time . 

The first record on each file* ccaisists of the year, number 
of records and number of files* in fixed decimal form. Each of the 
successive records contains the following information; 

Word 1: Initial time of record in hours frcm base time . 

(0.0^\jT December 5 I of year previous to launch) . 


* Pseudo file. Tape is prepared in overlapping two-year groups. 
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Equatorial coordinates of Ifercury in two-day intervals follow (9 x 
values, 9 y values, 9 z values) . Then 27 consecutive five-word 
blocks containing the equatorial coordinates, in four- day intervals of 


XVNE 

YVNE 

ZVNE 

XSE 

YSE 

ZSE 

XAS 

YAS 

ZAS 

XJS 

YJS 

ZJS 

XSAS 

YSAS 

ZSAS 

XUS 

YUS 

ZUS 

XNS 

YUS 

ZNS 

XPS 

YPS 

ZPS 

XBS 

YBS 

ZBS 


are followed by three 52-word blocks 
dinates of the 

XME YME ZME 


Venus with respect to the Sun 

Sun with respect to the Earth 

Hars with respect to the Sun 

Jupiter with respect to the Sun 

Saturn with respect to the Sun 

Uranus with respect t( the Sun 

Neptune with respect to the Sun 

Pluto with respect to the Sun 

Earth-Moon barycenter with 
respect to the Sun 

containing the equatorial coor- 
Moon with respect to Earth 


The Moon coordinates are stored in half-day intervals (O.O^, 12^0 UT) 
with distance measured in ER. All other tables are in AU. 

The equatorial coordinates of the planets and of the Moon are 
followed by their velocities, in exactly the same order. Moon velocities 
ai-e in ER/day. All other velocities are in AU/day. 

At present, an ephemeris tape is available for I965 - 19 ^ 9 ^ siiid 
1968 - 1982, written in 5 and 15 two year groups respectively, each of 
which overlaps one year. 
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C . Ephemeris in Core 

The astronomical tables are stored in core in 96 -bour intervals 
for the Sun eind the planets, and 12-hour intervals for the Moon. 
There are always 32-days of tables available, arranged in such a way 
that the value of time for which the interpolation takes place is 
not near either end of the table . 

In location TABLE(l), the time of the first entry from the 
initial time is stored. In TABLE(2) to TABLE(10) there are 9 
X coordinates of the Sun with respect to the Earth. The following 
chart indicates the storage locations of the remaining astroncanical 
data to be saved. 


TABLE(ll) 

to 

TABLE(19) 

y coordinates of the S’Jin with 
respect to the Earth 

TABLE(20) 

to 

TABLE(28) 

z coordinates of the Sun with 
respect to the Earth 

TABLE(29) 

to 

TA3IiE( 55 ) 

X, y, z coordinates of Jupiter 
with respect to the Sun 

table(56) 

to 

TABLE(82) 

X, y, z coordinates of Mars 
with respect to the Sun 

TABLE(83) 

to 

TABLE( 109 ) 

X, y, z coordinates of Venus 
with respect to the Sun 

TABLE( 110 ) 

to 

TABLE(136) 

X; y, z coordinates of Saturn 
with respect to the Sun 

TAEEiE(l37) 

to 

table(i63) 

X, y, z coordinates of Uranus 
with respect to the Sun 

TABLE(164) 

to 

TABLE(190) 

X, y, z coordinates of Neptune 
with respect to the Sun 

TABLE( 191 ) 

to 

TABLE(217) 

X, y, z coordinates of Pluto 
with respect to the Sun 

TA3LE(218) 

to 

TABLE(244) 

X, y, z coordinates of the Earth 
Moon Bary center with respect to 
the Sun 

TABLE( 245) 

to 

TABLE( 296 ) 

X, y, z coordinates of Jfercuiy 
with respect to the Sun 

table( 299 ) 

to 

table(363) 

X coordinates of the Ifoon with 
respect to the Earth 
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table( 364) to TABLE(428) 
TABLE(429) to TA3L(i.'93) 


y coordinates of the Moon with 
respect to the Earth 

z coordinates of the Moon with 
respect to the Earth 


Uiese are followed by the velocities: 


TABLE(494) 

to 

TABLE( 520 ) 

X, y, 
with 

z coordinates of the Sun 
respect to the Earth 

TABLE( 521 ) 

to 

TABLE(547) 

X, y, 
with 

z coordinates of Jupiter 
respect to the Sun 

TABLE( 5^) 

to 

table(6oi) 

X, y, 
with 

z coordinates of Venus 
respect to the Sun 

TABLE( 791 ) 

to 

TABLE(985) 

X, y, 
with 

z coordinates of the Moon 
respect to the Earth 


D. Perturbation Program 

The perturbation program solves three differential equations for 
XI, ETA, ZETA. The differential equation for XI, with the various 
terms replaced by the storages coniaining them, is representative of 
all three equations and is given below: 


D2XI * - GME 

[VCX)R( 1 ) /VCOR( 'O 

COMP( 1 ) /CX)MP( 4) ] 

- GMVN 

[VC0R(19)/VC0R(22) - 

C0MP(19) /C0MP(22)] 

- GMS 

[VC30R(13)/VC0R(16) - 

COMP(13)/COMP(i6)] 

- GMMR 

[VC0R(25)/VC0R(28) - 

C0MP(25)/CX)MP(28)] 

- GMJP 

[VC0R(31)/VC0R(54) - 

comp(3i)/comp(34)] 

- GMM 
+ OTHER 

[VC0R(7)/VCX)R(10) 

PERTURBATIONS 

COMP( 7) /C0MP( 10 ) ] 

where, for example, 

in the first term OE 

= is the mass o-^ ;,he Earth 


and VC0R(4) is the length cubed of the vector [VCOR(j), VC0R(2), VC0R(5)]« 
Fimilarly, in the other terms the denominator is the length cubed of the 
tor containing the corresponding numerator. In the case where the two 
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terms within each of the brackets are nearly eq.ual, they are computed 
by the special method described in Appendix E to avoid loss of accuracy. 

The contents of the COMP storage at any time, t, depends upon 
the reference origin at that time . 





CONTENTS 

OF COM? STORAGE 


Earth 

Moon 

Sun 

Venus 

Mars 

Jupiter 

Ref. 

Ref. 

Ref. 

Ref. 

Ref. 

Ref. 

C0MP(l) 

XVEO 

XME 

XSE 

X'.^ 

XMRE 

XJPE 

(l), (2), (3) = X, y, z 
(4), (5), (6) = R3,R,r2 

XEMi 

XVMD 

XSM 

XVMM 

XMRM 

XJPM 

(7), (8), (9) 

XES 

XMS 

XVSO 

XVNS 

XMRS 

XJPS 

(13), (14), (15) 

XEVN 

XMVH 

XSVR 

XWNO 

XMRVH 

XJPVN 

( 19 ), (20), (21) 

XEMR 

XMMR 

XSMR 

XVNMR 

XVMRO 

XJPMR 

(25), ( 26 ), (27) 

XEJP 

XMJP 

XSJP 

XVHKP 

XMRJP 

XVJPO 

(31), (32), (33) 

Here XVE 

refers 

to the 

X component of the 

vehicle 

with respect to the 


Earth, with corresponding definitions for the other quantities. An 
additional subscript of 0 denotes quantity derived frcmi the two-body 
problem. 


COKTENS OF VCOE STORAGE 

All 

Refs. VCOR(l) 

XVE (1), (2), (3), (4), (5), (6) = X, y, z, r3, R, R^ 

XVM (7), (8), (9) 

XVS (13), (14), (15) 

XWi (19), (20), (21) 

XVMR (25), (26), (27) 

XVJP (31), (52), (33) 
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XIII . Available Upon Request 


Numerous additions and improvenents are under development to the 
current OS/ 36 O version. These additions can be made available for 
general use to interested organizations A brief description of some 
of the additions are as follows; 

1. JLtiple Vehicles. 

It is now possible to integrate N trajectories simul- 
taneously (N = 1 to 30 ) . The user has the option of using the 
two-body solution for all trajectories or separate two-bodies 
for each trajectory. 

2 . Lambert ' s Problem . 

The program has the capability of generating its own 
initial conditions when one is interested in a specific inter- 
planetary trajectory. !Ehis option requires a starting Julian 
date, a desired flight time, and a target planet. Within this 
option, there is a further option which computes the initial 
conditions on the sphere of influence of the Earth or on a 
parking orbit inside the Earth's sphere of influence. 

3* J. P. L. Ephemeris. 

It is now possible to read the J.P.L. Ephemeris directly 
rather than by the method described in Section XI-1. This 
capability is obtained by adding a module of subroutines that 
would permit the trajectories to be integrated with respect to 
the mean equinox of 1950. 

4 . Small Corrective Thrust. 

5- Trajectory Search 

It is planned to automate the iteration scheme to go from 

**Address: Head, Theoretical Mechanics Branch, Laboratory for Space 

Physics, Code 643 ^ Goddard Space Plight Center, Greenbelt, 
Maryland 207?1. 
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two-body, to patched conic, to full trajectory, and to increase 
the number of variables to be adjusted, in optimal fashion. 
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Section XIV Methods of Integration 

The integration scheme employed by the Interplanetary Trajectory 
Encke Method program is a sixth order backward difference scheme, 
initiated by a Runge-Kutta scheme. The routine used is Nevton- 
Gregory integration scheme for general second order difference equations. 
( See Appendix C . ) . 

The program has an option (N0PT(71) ^ O) for integrating Beta 
instead of time. (See Appendix M.) . Computation is much faster in 
this mode, however, the user is cautioned to choose delta beta with 
caro . 


The program has an option (N0PT(20) =» 2 or J) for integrating 
second order differential equations by means of intemolating a table 
of second order derivatives . The size of the vable and the option of 
using predict or- correct or, or predictor only, are Inputs to the program. 
(See Appendix j) . 
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MATHEM/ITICAL APPENDIX 


A. INTROmCTION 

The problem of orbit determination over long time periods 
requires a precise technique for integrating the equations of 
motion. Reference 4 contains an analysis of em integration 
procedure that yields, the mlnimim loss of information due to the 
accumulation of numerical round-off errors. The Encke perturba- 
tion method has been shown to require minimum machine computation 
time for a minimun loss of numerical accuracy. The orbit pre- 
diction scheme presented herein uses a modified form of the Encke 
method with the initial position and velocity vectors replacing the 
conventicmal P and Q vectors of the Encke scheme. ^ avoiding 
reference to the position of perigee, it is possible to avoid 
numerical ambiguities arising from near-circular orbits and orbits 
of low inclination. 
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B. EQUATIOIE OF MOTION 

In a Newtonian system, the equations of motion of a particle in 
the gravitational field of n attracting bodies and subject to other 
perturbing accelerations such as thrust, drag, oblateness, radiation 
pressure, etc. are given by 



(B.l) 


These equations are put into observable form by referring them to a 
reference body c. The equations of motion of the reference body are 


•• 



n 


V 

L 


i=l 

ij^c 



(B.2) 


Subtraction of Equation (b. 2) from Equation (b.I) results in the 
equations of motion of the vehicle with respect to the reference body c. 







r 


VC 







J 


(B.p) 
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C. IMTBGRATION PROCEEUEE 

If Equation (b.J) is integrated directly by sane numerical 
scheme, there results, after a number of step-by-step inoegrations 
an accumulation of error which leads to inaccurate results. To 
avoid this loss in precision, it is convenient to write Equation 
(B.3) in the form 


•# •• 

The velocity and displacement vectors can be written as 

■>vc ’ \ ^ 4E 


The reference body (the one in whose sphere of influence the 
vehicle travels) is chosen so as to minimize the perturbations. 



Equations (C.4) constitute the equations of motion of the Kepler 
T'l'oblem and are solved as described in Appendix D. 


(C.l) 

(C.2) 

(C.5) 


(C.lf) 

(C.5) 
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Equations (C. 5 ) are integrated numerically. The integration 
scheme employed by the UTM program is a sixth order backward difference 
scheme , initiated by a Runge-Kutta scheme . The routine used is a Newton- 
Gregory integration scheme for second order difference equations written 
by S. Pines and J. Mohan of Analytical Jfechanics Associates, Inc. 

As derived in Appendix D, the solution of the Kepler problem may be 
represented by the vectors Rq, Ro; the scalar a and the rectification 
time to* 

The rectification process consists of moving R , R into the 
. vc' VC 

locations Rq and t into to and the computation of a and n. 

For conputational convenience, the coefficients appearing in 
Equations (D.2) are also conputed during rectification. 
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D. SOLUTION OF THE KEPLER TWO-BODI PROBLEM 


The unified formulation of the two-body problem is used for 
both elliptic and hyperbolic cases . 

e ’,/N • ® 


» =■ (a) 

Fa(of) • I - fij + 




00 


I 


I 


(-0) 

(2i + 3 ) I 


( ~or) 

(2i + 2) I 


i»0 


(D.l) 


FaCor) =» 1 - or Fj^ 
F4(o) = 1 - o F2 


and 

f =» 1 - — F 2 
ro 


g 


— P Fa 



f ’ “ ^ P Fa 
g * 1 - p2 F 2 


(D.2) 
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where 


do = Ro * Ro 


r = F 2 + ro F 4 + ~ 3 Fa 


/m- 


R = f R + g R 
o o 


R = f Ro + g Ro 
o 


2 

V_ .-I 


(-■ — ) 

\ ro |i / 


Vo = Ro • Ro 


a is determined from the modified Kepler equation 


_ 

v5T At = Fj + ro 3 F 3 + — F 2 


(D. 3 ) 


See Figure 1 for the two-hody orbit which results fran the solution of 
Equation (C.4) with the initial conditions: 


Rj^(to) - Ryc^'to) » Ro 

Rj^(to) ■ ^ 


(D.if) 
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y Axis 



Figure 1. Geometry of the Elliptic Two-Body Orbit 
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E. COMPUTATION OF PERTURBATION TERMS 


The terms accounting for the Encke term and the planetary per- 
turbations appearing on the right hand side f Equation (C.5) 

P R 

involve numerous terms of the form -5 - where R and R© may 

r 1*0 

differ only by small amounts. For the Encke term, for Instance 
R - Ro = § which is small, and for the planetary perturbations, the 
difference is R^^ which also often is small. 

A computation scheme, which avoids loss of precision due to 
the subtraction of nearly equal terms and which also is correct 
when R is not small, is en^loyed. This scheme is described 
below; Find 



u 



( Rq *1 ar) * AR 


R 

“3 

r 




R(u® 

T 


+ 3^^ 5u) 


1 + 


) 


o 


(E.l) 
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CONCLUSIONS 


The method presented yields accurate trajectories using 
relatively little computer time . Summarizing some of the important 
features : 

1 . All significant solar system bodies may be included Wi.thout 
undue complications . 

2 . Since the perturbations only ar^ integrated, the allowable 
integration interval is fairly io-rge over most of the path. 

Even in the vicinity of Earth or another planet a relatively 
large interval ( compared to other schemes) nay be used 
without limiting the stability and accuracy of the solutions. 

3 . The perturbations are kept small in two ways. First, the two- 
body orbit is rectified whenever the perturbations exceed a 
specified maximum value con^jared to the corresponding unper- 
turbed values. This limits error build-up with respect to 
particular reference body. Second, the reference body of the 
two-body problem is changed from Earth, to Sun, to planet 
accordingly, as that reference body would contribute the 
largest perturbing force otherwise . 

4. This method will handle circular orbits, zero inclination, etc. 
The problem is defined in terms of parameters which have real 
physical significance (namely, the position and velocity 
vectors) which are directly relatable to measurable quantities . 



OBLATENESS TERMS 


A subroutine called GOBL was modifed to obtain, inertial 
Cartesian coordinates, the gravitational perturbation acceleration 
due to a rotating nonspherical body whose mass coefficients are 
given in terms of the zonal and tesseral harmonics. The method 
described herein avoids the classical singularity, which occurs 
for polar passage, when using spherical coordinates to describe 
the gravitational potential. G?hls method also minimizes the 
numerical error incurred on double appllcetion of coordinate 
transformation from inertial to body-fixed and back again in the 
case of a rotating nonspherical gravitational body. 



DERIVATION OF THE EQUATIONS 


The potential at a point R, in the coordiniate system fixed in the 
body, is given by 

CO 

“P * • I l,f) - I * ®n,m ■"'>]} 


n=»2 


msl 


where 


z 

u = — 
r 


tan X 


1 

X 


Pn(u) 


2°n: du 


d / 2 

(u 

n 


- 1 ) 


( 2 ) 


P (u) = (1 - — P (u) 

m Wl ' ' ' ' «M ' 


n,m 


du‘ 


m n 


The accelerating force vector, in the body- fixed Cartesian coordinates, 
is given by 


where 


F = ^ ^ ^ + 1? lii- + is iL. 

9x^ Br dx^ dx^ 9X dx^ 




R 




R 


dX 

dx^ 


1 

r( 1-u^) 


k X R 


k = 


0 

0 

1 


( 5 ) 


ih) 
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Combining the scalar coefficients of the vectors R, k, and k x R, we 
have 

™ .. h 


F = 




n=2 


n n 


^ ^ r ^ 

- ) — i ) [u P’ +{n+l)P ][C cos mX + S sin mx] 

L n+2 i L. n,m n,m-'‘- n,m n,m . 


n»2 


m=l 


j 


|Jk ^ 

- y — P' - y P* (C cos mx + S sin mx)l k 
L. _n+2L n n L. n,m n,m n,m J 


n 


( 5 ) 


n=2 




n n 


+ ) > P (S cos mX - C sin mX) kxR 

Lj ^^2 ^ (l-u^) ■* 

n»2 r m=l ' " ' 


Examination of Eq. (5) indicates that, for a polar passage with u « 1, 

A A 

a singularity occurs in the coefficient of the kxR term and in the 

derivatives of the associated Legendre polynomials P (u) . In order to 

m 

remove this singularity euad to avoid numerical inaccuracy in trajectories 
close to u = 1, it is convenient to change the coordinates from the 
spherical r, u, X system to the four-parameter system defined below. 


Let the body-fixed Cartesian vector be 

s 

R = r It 
u 


( 6 ) 



where 


X 

s = — 
r 


t 

r 


(6a) 


z 

u =- 


In this coordinate system, the potential, cp, is given by 


00 


n 


■P = r - I 0) ^ = . S„^„l(s,t)]}] ( 7 ) 


n=*2 




where 


n+m 

« id ^2 

A =» — ( u 

n,m ^n , , n+m 

' 2 nl du 


- 1 )“ 


R (s,t) * Real Part of (s + i t)°^ 
m 


(7a) 


Iuj(s,t) = Imaginary Part of (s + i t)°^ 


i 


In this system, the acceleration force vector in body-fixed 
Cartesian coordinates is given by 

T? = ^ ^ M ^ ^ ^ ^ 

9r Bx^ ^ ds dx^ 3t 3x^ Bu Bx^ 
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where 


■ r - ® 


3s ^ 


1 ^ 
— X 

r 


a~t 1 ^ 

a — i 

3x^ r 


r 




(8a) 


3u 

3x. 


ii; 

r 


r 



’i 

A 

d 

A 

► « 
0 

m 

0 

0 

9 

; J = 

0 

L J 

; k = 

0 

1 

L J 


Combining the scalnr coefficients of the vectors R, i, j suad k, we 


have 


F,(2s!.l^.i^.“afriR+i^i + i^3+i22£ 
\ar r 3s r at r ^u/ r as r at r au 


(9) 


or 


F = or^R + a^i + cr^O + 


(9a) 


For the perturbation coefficient of R, we have 


|J^ £L 

y — ^(jPaA -+(n+l)A 1 - 7 fu A ^^-Kn+l)A 1(R C +I S ) 

L n+2 L nL n,! n,oJ L L n,m+l n,mJ m n,m m n,m 


n=2 r 


m«l 


‘ y “ ^ (s I T+ t R^ -)S„ 3 

Li n^QiL m~l m-1 n^m m-1 m**! n^oLU 


m»l 


(lo; 
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However 


s R 

. - t 

I , 

= R 

m-1 

m-1 

m 

s I 

, + t 

R 1 

=» I 

m-i 

m-1 

m 

two inner 

summations 


to lead 


n 


' y [j + (n+m+l)A „ 1 (R C „ + I ) 

Li L H;ia+X n,m J m n^m ni n^m 


m=l 

Furthermore, if we let 


J = - C 
n n,o 


the expression for the coefficient of R may be written as 
_ n n 


O' = 
r 


= ) fuA +(n+m+l)A 1 (R C +IS ) 

L, n+2 L, L n,m+l n,m J m n,m m n,m 


n=2 m*0 


A recursion equation for A^ ^ may be derived, which yields 


^ \ ^ 4-1 (n+m+l)A 
n+l,m+l n,m+l n,m 


The final compact expression for the coefficient of R is given by 


n n 


“ [J, a “ 

“r= J 


n=2 


ir=0 


( 11 ) 


(lOa) 


(10b) 


( 12 ) 


( 13 ) 
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Turning to the coefficient of k, we have 


« a 


n 


03 


= - y — -T fj A . - y A (R C + I S ) 
Zj n+2 L n n,l u n,m+l m n,m m n,m J 


n=2 


m=l 


Using the convention that J = -C , Eq. (I 5 ) becomes 


n n,o" 


00 0“ 

|i a 


as 




n=2 m=0 


The coefficients of 1 and 3> respectively, follow without modification. 


n n 




n=2 m=l 


and 


n n 


QT2 


’ y y “ n. 1 

Z_i n-r<j n,m L ni"l n,m m-1 n,m J 


n=*2 n»l 


( 14 ) 


(l4a) 


(15) 


(16) 
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RECUBSION EQUATIONS FOR THE M3DIFIED LEGENDRE POLYNOMIAI£ AND THEIR 
DERIVATIVES 


We seek a recursion equation for the modified Legendre polynomial, 

A . We have 
n,m 




n,m 


( u) = — - P^( u) 


du 


m n 


(17) 


In any standard reference on Legendre polynomials, we may obtain the two 
recursion equations. 


(n+l)P + u P' = P' _ 
n n n+1 


( 18 ) 


and 


(2n+l)P + P ' = P’ , 

n n-1 n+1 


(19) 


In terms of A (u), Equations (18) and (19) become 


(n+l) A +uA - =A_- 
n,0 n,l n+1,1 


(18a 


and 


(2n+l) A +A 
' n,0 n-1,1 n+1,1 


(19a. 


Combining Eqs. (l8a) and (l9a) and eliminating A , we solve for 

3^ ^ U 


An+i ^(u) as follows; 


^ (2+i)P* - (1+^)P' T 
n+1,1 n+x n n n n-1 


T> • 


( 20 ) 
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or 


A - u(2 +^)A - (1 +^)A^ , ^ 

n n j> i!i 

Successive differentiation of Eq. (20a) yields 

A = (m-1) (2 + -i- V - -,•*■ u(2 +~ )A , -(1 + -~ V o 

n,m \ n-1 / n-l,m-l \ n-1 /^n-l,m \ n-1 /^n-2,m 


This is the required recursion equation. Furthermore, we have 


A ^ 

0; 

m > n 

n,m 



II 

a 

1 • 

5.5.... (2n-l) 

A , 

=! U 

A 

n,n-l 


n,n 


Starting with Eq . ( 22) , we may generate each row of A , for fixed 

n ^ in. 

n, by retaining the two previous rows of A^ ^ and through application of 
Eq. (21) for m = n-2, n-^, , , , to m=l. Thus, only three rows of A^ ^ 
need be retained for a given a- 

To obtain the result given in Eq. (l2), we need only differentiate 
Eq. (18a) and obtain 


n+l,m+l 


u A 


n,m+l 


+ (n+m+l)A 

n ,m 


An alternate recursion formula 
(21); is 


for A , which is more stable than Eq. 
n,m^ ^ 


n.m 


n - m 


(u A 


n,m+l 


■^n-l,ra+l^ 
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(20a) 


( 21 ) 


(22) 


( 12 ) 


(21a) 



RECURSION EQUATIONS FOR R AND I , AND THE ACCELERATION EQUATIONS FOR 

m m 

A ROTATING BODY 


For the zero power of s + it, we have 

(s + it)° =* 1 


(23) 


Thus 

Ro * 1 

lo * 0 


(23a) 


From Eq. (ll), we obtain the recursicai equations for R^ and I^. Given 
s and t, we have 


We may define 


Rjs.t) 


s R - - t I T 
m-1 m-1 


s Vi " Vi 


D 

n,m 

E 

n,m 

F 

n,m 


R (s,t)E + I (s,t)S 
m n,m m ' n,m 


R^ ,(s,t)E^ „ + 1 i(s,t)S^ ^ 
m-1 n,m m-1 n,m 


R ,(s,t)S - I ,(s,t)C 

m-1' ^ n,m m-1' ^ n, 


m 


(11) 


( 24 ) 


The final desired form of the equations for and 

are given by 
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(25) 


In the body- fixed system, we have, for the acceleration. 


F =* cy^ R + i + (>2 «) ^ 


(9a) 


Let the rotation matrix N(5 x 3) represent the transformation of a 
rotating body from inertial to body- fixed coordinates . 


R =» N R 


inert 


( 26 ) 


Then the inertial acceleration is 


T 

F » N F 

^ inert 


(27) 


Or, in the inertial system, we have 


F . ^ - a Rj . + Of-. N or-. N., 

inert r inert 1 2 2 3 3 


(28) 
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H. TRANSFORMATION EQUATIONS FROM GEODETIC POLAR COORDIMATES TO 
CARTESIAN COORDINATES* 

The geodetic polar coordinates in the program ai« referred to 
an ellipsoid of revolution. The equation of a cross section is 
given by 


where 



1 


=* a^( 1 - e ’) 


The slope of the normal, along which h is measured is given by 


1 Si^Z 

tan cp ( See Figure 2) 

^ b^x 
dx 


(H.l) 


(H.2) 


and 


z Id" ^ 

tan cp' » — « tan cp » ( 1 - e^) tan cp 

b^ 

Eliminating x between equations (h.I) and (H.2) and solving for 
z results in: 


a(l - e^) sin cp 
(l - e^ sin^ cp) 


*For geocentric (i.e. e^ « O) polar coordinates, c » s = 1. In this 
case the latitude input is interpreted as declination. 
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Figure 2. Relation Between Declination, Geocentric and 
Geodetic Latitudes 


and from eq.uation (H.2) then 


a cos (p 

(1 - sin^ cp)® 


H-2 


In units of a. R and R are then given by equation (H.5) 
e « (1 - e^ sin^ cp)”^ 


s = (1 - e^) c 

X - (c + h) cos 9 cos (0 - 0o) 

y = (c + h) cos cp sin (0 - 9 o) 

z = (s + h) sin cp 

cos (0 - 0 o) 

(H.3) 


i u v{(slnY COS Cp - COS Y cos A sin cp) 
“ cos Y sin A sin (0 - 80 ) j" 


V “^(sin Y cos cp - cos y cos A sin cp) sin (0 - 0o) 


+ cos Y sin A cos (0 - Qq)]" 


V sin Y sin 9 + cos y cos A cos cpj 


These equations include the effect of the rotatican of the earth. 

The longitude of the vernal equinox ( 0 ©) at launch tinse is conq>uted 
by the program from Newcomb's formula. 
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I. TRANSFORMATION EQUATIONS FOR RADAR SIMULATION 

The program computes sight angles (in an azimuth-elevation 
system) , slant range and range rate data for up to 30 radar 
stations . Die vehicle coordinates are transformed from a system 
of geocentric cartesian coordinates (xyz), the x-axis in the 
direction of the vernal equinox and the x-y plane in the equa- 
torial plane of the earth to the required topocentric azimuth 
elevation system. This is accomplished by a series of coordinate 
transformations as follows: 

1 . A rotation of the coordinate system about the z-axis 
through an angle RA^ so that x y plane is in the meridian plane 
of the station. 


x' 

= X 

cos RA 

s 

+ y sin RA 

s 


y’ 

=» -X 

sin RA 

s 

+ y cos RA 

s 

(I.l) 

z' 

K Z 





The velocity transformation must tsike the rotational velocity of 
the new coordinate system into account . 

x' » y' (B + X cos RA + y sin RA 
•' e s s 

y* = -X* (jD - X sin RA + y cos RA 
6 s s 

z* = z 

where x*, y', z' are the rotated coordinates and RA is the right 

s 

ascension of the station and od is the sidereal rate of the earth's 

e 

rotation. The G.H.A. necessary to obtain RA^ from the station 
IcMigitude is con5)uted by the program. 
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2. A translation of the origin of the coordinate system from the 
center of the earth to the station in question 

x'* = x* - (c +h) cos cp 


z'' = z' - (s +h) sin cp ( I *2) 


where 

c ^ (l - e^ sin^ cp) ^ 
s = (l - e^) c 

x'*«x'; y''=y'> z''»z’ 

where x' y' S z' ' are the translated coordinates, cp is the 
geodetic latitude and h the height above sea level of the station 
in question. 

5- A rotation of ( 90 - cp) about the y' ' axis to place the (x*',z*') 
plane into the horizon plane 


X I « I ^ X ' ' sin cp + z ' ' cos cp 

yl f I = y' « 

z 1 1 • _ x' ' cos cp + z ' ’ sin cp 

X ' ' ’ = X ' ' sin q) + z ' ' cos cp 
y' ' ' => y' ' 

z ' ' ' 3 - X cos cp + z ' ’ sin cp 


( 1 . 3 ) 
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Now x''', y' ' ' , z'*' are the coordinates of the vehicle in a topo- 
centric azimuth elevation system, with z ' * ' axis pointing to zenith 
and the x’*' pointing south along the meridian. Range, range rate, 
azimuth and el-' vat ion are then given by 


= (x' 


I t i2 


)i = 


Slant range 


( 1 . 4 ) 




Elevation 


( 1 . 6 ) 



A 


r tt - A* A' < TT I 

\ 5tT - A* A* > TT J 


(1.7) 


1-3 



J. A VARIABLE OREER INTERPOLATION SCHEME 

FOR INTEGRATING SECOND ORDER DIFFERENTIAL EQUATIONS 

A subroutine called AMAINT was programmed to integrate second 
order differential equations by means of interpolating a table of second 
order derivatives. The size of the table and the option of using 
predictor-corrector, or predictor only, are inputs to the program. This 
subroutine operates in a fixed step- size mode only. 

The program uses a self-starting scheme instead of the usual 
Runge-Kutta starter to build a table of second derivatives. This scheme 
employs the following technique; a first guess for the tables is made by 
stepping up the independent variable and calling the derivative routine; 
now, using this table and the predictor only formulas, calculate every 
point on the table in succession. When either a) all of the first 
derivatives, or b) all of the second derivatives, from one iteration to 
the next agree to 15 digits, we consider the scheme converged and we now 
have a starting table. Normally this converges in 8 iterations. For a 
twelfth-order integrator, our starter routine calls the derivative 
routine 89 times; the Runge-Kutta starter would require 176 calls. 
Whenever a step size other than the normal is required, the program can 
take that step using the stored table of second derivatives, rather 
than using Runge-Kutta again as is usually done. 
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Derivation of Equations Used 

The function to be integrated is: 

y(t) = f[y(t), y(t),t] (l.o) 

We assume that we have n function values 

^s-n' ' " ^s-1 

over some constant increment h of t 

I — , 

h 


t 


s-n 


^s-n+1 


t 


s-1 


+ 



where s is chosen to be at the midpoint or close to it . We choose s 
in such a manner that our coefficients will be integers for a large n. 
Q5ie maximum n depends on the word size of the caaputer being used. 

The function may be extrapolated to t using an nth order Lagrange 

s 

interpolation formula : 


s-1 

n (t - t.) 

s-1 i = s-n 

LjLJ 




J=s-n 


n 

i =» s-n 

i J 




( 1 . 1 ) 
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This results in a pol^Tiomial in t 


y(t) 


n 0 
/ 'T' 

=» / ( '/ S , 

Li \ j 

j5»l i=:n-l 


,1 ' 


S. = 1 
J,n-1 


Integrating we get; 


y(t) 


” ° s. 


i=n-^l 


n 


0 e 


= Z [( I (iii)(n- gT * =j ^ 

Jal i=n-l 




^ r ^ s s-i> 

L L (i+lj (i+2) 

j=l i=n-l 


(t - t - 
s s-1 


° s. 




i=n-l 




I[ I 

ji*! i^n-1 


fl+1)" /j ^ 


Consider each interval as 1 since they are equally spaced. Now 

t=h’S. S. .isan nxn matrix S*D where S is formed in the fol- 
s j,i nxn nxn 

lowing manner; The first row is all one’s; the second is (-) the sum 
of the t^'s of Equation (l.l) in the numerators of the varying j's; the 
next is the sum of the products of the t^'s taken two at a time; etc. 
Thus we have a matrix of the coefficients of the polynomials formed by 
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( 1 . 2 ) 


(1.5) 



the numerators of (l.l) sitting column-wise . The D matrix is a 
diagonal matrix whose elements are the denominators of (l.l) . 


We now form matrices A and B where 


A = a_, , 
nxn ij 


( s-n-l+i+e, ) 
k 


n+l-J 


- (s-n-l-fi)^'^^’’^ 


n+l-j 


and 


B = . 

nxn ij 


(s-n-l+i+ej^)“'^"‘^- (s-n-l+i)^'^''*^- (n+2-j)ej^( s-n-l+i) 


Now 


ASD =» or. 




BSD = 


Equation (l.2) becomes 


y(t^) = (y(t^.^) y(t^.j^) + ^ 


J*1 


and Equation (I. 5 ) becomes 


n 


y(t^) » y(t^.^) + h I 

J»1 

where e = li,/h, h^ ^ h when t -t h. These are the predictor 

K «L X S S*X 

equations . 

If we desire a corrector formula, we generate 
c c 

. and B . J » 1, n+1 

■“yj n,J 

We use Eqs. (2.l) and (2.2) to find y(t ) and y(t ) . Then using Eq. 


J-4 


n+l-J 


( 2 . 1 ) 


( 2 . 2 ) 



(l.O), we find y(t ) . Now we have a table of n+1 f 

® 0 
our corrected functions and first derivatives by: 


s and we can find 


n+1 


I Pn,j ^ 


j=i 


n+1 


yc'*s> = * •' Z “n,J 


J=1 


When h^ / h, the corrector formulas are not valid. Qherefore, when an 
odd integration step is required, predictor only may be used. 


The and matrices are used to generate the initial table of 

f ’s at each restart by means of an iterative method. In generating 
^ c c 

the coefficients P . on a computer, we attempt to keep 

then integers by limiting n to conform to the word si?.e of the computer 
being used. We also use a least ccsnmon denominator technique to perform 
one, instead of n, divisions. Thus we see that the word size of the 
ccmputer controls the order of the integrator to be used. When our 
constants do not sit as integers in the machine, the resulting round-off 
causes biased errors. 


A study was made using this routine on an IBl 36 O Model 9 I computer 
to integrate the sine and cosine functions. The Uni‘'mc UO 8 computer 
was also used to runa a small number of cases. Orders from T to 12, and 
integration intervals from tt/i6 to tt/526, predictor only and precitor- 
corrector, were used with the following results: 
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n - Number of points (not differences) 


Predi ct or- Cor re ct or 


Sin TT 

n = 6 

tt /32 

.68 X 10“^ 


tt /64 

.27 X 10"“ 


tt/128 

.12 X 10"^^ 


tt/256 

.28 X 


tt/512 

.55 X 10"^'^ 

n = 7 

rr/32 

.75 X 10"^ 0 


tt /64 

.40 X lO"^^ 


tt/128 

.12 X 10“^ 3 


17/256 

.29 X lO"^'^ 


tt /512 

.54 X 10 ‘‘ ^ 

n = 8 

n/52 

.58 X 10"^^ 


n /64 

.49 X 10 "^“^ 


tt/128 

.16 X 10"^ 


tt/256 

.26 X 10’^“^ 


tt /512 

.49 X lO"^'^ 

n = 9 

tt /52 

.12 X lO"^^ 


tt /64 

.18 X lO"^*^ 


tt/128 

.17 X lO"^*^ 


n/256 

.50 X 10 "l^ 

n ^ 10 

tt/52 

.52 X 10"^ 2 


tt /64 

.11 X 10"^ 


tt/128 

.19 X 10'^^ 


tt /256 

.51 X 10"^ ^ 
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Pr e di c t or- C orrector 




Sin TT 


n ^ 11 

tt/i6 

.12 X 10“^^ 


tt/52 

.84 X lO"^^ 


tt/64 

.11 X lO"^^ 


tt/128 

.18 X lO"^^ 

n =» 12 

it/i6 

. -11 
.54 X 10 ^ 


tt/52 

.81 X 10‘^^ 


tt/64 

.74 X 10“^^ 


tt/128 

.71 X 10"^^ 

Predictor Only 



n =* 6 

tt/52 

.42 X lO"^ 


tt/64 

.69 X 10"® 


tt/128 

.11 X 10‘^ 


tt/256 

.17 X 10”^^ 


tt/512 

.21 X 10"^^ 

n = 7 

tt/52 

.84 X lO"® 


tt/64 

.45 X 10"^^ 


tt/128 

-12 

.19 X 10 


tt/256 

.18 X 10"^^ 


tt/512 

.45 X 10"^^ 

CO 

It 

a 

tt/52 

.57 X 10"^° 


tt/64 

.15 X 10"^° 


tt/128 

.62 X lO"^^ 


tt/256 

.51 X lO"^^ 


tt/512 

.51 X lO"^^ 
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Predictor Only 

At 

Sin rr 

n =* 9 

tt /32 

.82 X 10“^° 


tt/64 

-12 

.11 X 10 


tt /128 

.U X lo’^^ 


tt/256 

.25 X 10“^^ 


tt/512 

.46 X lo’^^ 

n = 10 

tt/i6 

.16 X lO"^ 


n/32 

,31 X 10"^° 


Tx/6k. 

.34 X 10"^^ 


CO 

.15 X 10'^^ 

n = 11 

tt/16 

.48 X 10'^ 


tt/32 

-IP 

.73 X 10 


tt/64 

.66 X 10 


tt /128 

.59 X 10'^^ 

n « 12 

VO 

.53 X 10"^ 


T" ' 

.23 X 10'^ 


tt/64 

.40 X 10"^^ 


tt/128 

.41 X 10*^^ 
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Uni vac 1106 OiJy 


Predictor-Correct or 


Sin n 

n = 12 

tt/i6 

.34 X 10“^ 


tt/32 

.17 X 10“^^ 


n/64- 

.27 X 10“^^ 


tt/128 

.29 X 10“^^ 


From the above we see that, as we increase n, we can also increase 
the At up to a point . The optimum appears at n 3 u for the XBM ^60 
Model 91. On the Unlvac U08^ the optimum was at n - 12. When using ;he 
predictor-corrector method^ we can maintain the same accuracy as pre- 
dictor ooly^ using twice the integration interval at the c^imum n. 
Qlierefore, we conclude that using the predictor-corrector judiciously 
(i.e., with the proper integration interval for the function to be 
integrated) may prove to be better than predictor only, with a very 
slight cost in computer time. 

Our optimum At for the sine function was tt /64 for predictor only 
and tt/52 for predictor-corrector. Bie preceding data ar^ true for the 
ISM Model 91 cooQuter suid for the functicHi tested. Ve believe that, 
with a computer having a larger word size capacity, we would find a 
larger n to be optimal since it would be possible to generate Integers 
Tor the integration coefficients. 



A change in the hardware of the lEM 36 O Model 9 I was made which in- 
creased precision in multiplication. Using this machine, we ' len ran 
trajectories similar to those discussed in the Fehlherg paper. ^ Instead 
of using a rotating coordinate system, we used an inertial coordinate 
system and integrated in an Encke mode, using 3 as an independent 
variable instead of time. We converted Fehlberg's initial conditions 
(xo = 1.2, io = 0, yo =» 0, yo =» 1.04-9557509850520, = 1 / 82 . 45 ) to our 

system and they became 

xo = 1.212128562765511 

xo = 0. 

yo =* 0. 

yo • .162771052954^991 

. 1/8S .45 

Hg = .987871447254689 

(03ie o subscript denotes the Moon's initial conditions.) 
m 

The starting point of our trajectories was at apogee behind the Moon, 
on the line of centers of the Barth and Moon. Our integrations traced a 
figure eight around the Moon and terminated at the second apogee . We 
ran six different trajectories, using 5 different AP's with 11th and 12th 
order predict or- corrector sehenKs for each. The results were as follows; 


1 "New One-Step Integration Methods of High-Order Accuracy Applied to 
Some Problems of Celestial Mechanics", Erwin Fehlherg, NASA George C. 
Mir shall Space Plight Center, Huntsville, A1b» 
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A 3 No. 

of Steps 

Order of 
Integration 

Starting Apogee 

2 nd Apogee 

.0078125 

1258 

12 

1.2128562765311 

1 .2128562765114 

.0078125 

1258 

11 

1.2128562765311 

1.2128562764704 

.00590625 

2515 

12 

1.2128562765311 

1 .2128562765340 

.00390625 

2515 

11 

1.2128562765311 

1.2128562765550 

.001953125 

5028 

12 

1.2128562765311 

1.2128562765353 

.001953125 

5028 

11 

1.2128562765311 

1.2128562765368 


For this particular problem, A 3 = .00390625 with the 12 th order 
predictor-corrector scheme appears to be the optimum trajectory. 
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Use of Routine: 

In order to utilize this subroutine, a programmer must supply the 
following common cards: 

cx)m«)n/ami/i3eilt, t, nn, xi(50), xid( 50), 12x1(30), ifst( 5 o). 

COMMDN/AME/m( 14) , IS, Nl, N2, N21, NBQN. 

■j?h' . non names may not be changed. 

The variable names may be changed providing the notation for integer or 
real numbers is maintained. 

M(14) is a dummy block of storage used internally by the routine. 

IS is the reference point for generating the integration coefficients. 
K1 is the order of integration. 

^ Nl+1 is used for the corrector. 

N21 =* 0 is used for the predictor only. 

N21 =1 is used for the predict or- corrector. 

NBQN is the number of equations to be integrated. 

DELT is the normal integration interval. 

T is the independent variable . 

DTI is the current integration interval, oddball or normal. 

Xl(l) is the i^J^ function value. 

XID(i) is the i"th first derivative. 

D2Xl( l) is the i'*'^ second derivative . 

IFST(I) = 0 means the ith equation is second order. 

IFST(I) = 1 means the i"^^ equatiaa is of first order. 

For COMMON AMI, the user must define Nl, N21, and NBQN in his program. 

For COMMCN AML, the user also muot define DELT, DTI, Xl(l), XID(l), 
D2XI(I), IFST(I) (I 1, NEON) at T = To* 
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The user must also supply a subroutine called DERIV, with the 
proper COMMON included and which computes Dexi( l) as a function of Xl( l) , 
XID(i) and T. If first order equations are desired, the first deriva- 
tive should be stored in Dexi(l), and Dexi(l) should be a function of 
XI(I) and T. 


The following calls to AMAUn? are usel: 


Call AMAINT ( O) 
Call AMAINT ( 3 ) 

Call AMAINT (l) 
Call AMAINT (2) 


initializes the integrator and calls LERIV. 

generates integration table. OET.T may be changed 
before this call and the user must make sure that the 
proper Call DERIV has been made. 

takes one normal integration step and calls DERIV. 

requires th' 0 DTI be defined prior to the call, takes 
an integration step of size DTI, and calls the 
derivative routine. 


J-13 



K. DRAG COMPUTATION 

The drag acceleration is computed, assuming a spherically symmetric 
atmosphere rotating with the earth. Thus: 


D = - 1/2 p A 



(K.l) 


where 

Veff =* R - u) X R 


u) is the sideral rotation rate vector of the Earth. 
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L. COMPUTATION OF SUBSATETJ.ITE POINT 

The geodetic coordinates of the subsatellite point are ccaaputed 
by the following method: 

The geocentric latitude (declination) is obtained from 


sin cp* 


z 

r 


This latitude is then corrected ^ geodetic latitude by the formula 

cp - cp’ + a .2 sin 2cp' + a 4 sin 4.cp' + ae sin 6c^- ag sin 8cp’ 

where 



(L.l) 


(L.2) 


(L.3) 
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e = the ercentricity of the earth 
r w the distance from earth's center 
See Reference 5* 

The geodetic height is then given by 
h = r cos (cp - cp‘) - y l-e^sin^9 


(L.4) 


The longitude is obtained by subtracting the sidereal time of Greenwich 
from the right ascension given by 

tan RA = (L.5) 


L-2 



M. CHANGE OF INDEPEljriEIIT VAEIABLE - BETA MODE 


According to the standard Encke method, we introduce a differential 
equation 



(Y.l) 


In the construction of the closed- form solution for (Y.l), a parameter 
0 arises. It is related to t by Kepler's equation, 

t = to (Y.2) 

/m- 

where f is a transcendental function of 0 and is obtained by summing 
several power series. 

If t is taken as the independent variable. Equation (Y.2) has to 
be solved for 0 by an iterative method, requiring numerous time- 
consuming evaluations of the function f for each integration step. 

Using 0 as the independent variable, however, only requires a single 
evaluation. 


It remains, of course, to see what becomes of equation (Y.l) and 


^ = X - p 






(Y.5) 
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If p is the Independent variable. We have, from Kepler’s equation, 
that 


« = id 


at any point along the solution of (y.I) . Thus 


(Y.4) 


P 




at any point along the solution of (Y.I) and the initial conditions 
become 

Xo|xo| 

p(Po) = xo and p'CPo) = — - - 


when 


00 0(to) = 0 

Now the solution for (Y.I), p and p*, can be written in closed form for 
any p. As auxiliary quantities in this solution, we have |p| and 

D * P ’ 

D = . They are ccmaputed as functions of p before p and p ’ are 

Ipi 

known; that is, with accuracy at least as good as that of p and p' . 

Not only are they needed and easy to compute, but they also have the 
interesting property that 


dt _ lol 

(Y.5) 

d^ ^ __D_ 

Mr 2 


and 



Thus Equaticn (Y.l) is solved more economically in terms of p than in 
terms of t . 


Now we turn to Equation (Y« 3 ) • To treat it, we want to express 
' in terms of From (Y» 5 ) have that 


I' 




Differentiating with respect to p. 


I 


I t 




(Y.6) 


•J M. 


+ ? 


, G. J. 


- IpI^ ( 


).m! 


F + 

Ip 


(y.7) 


Thus (Y.7) is the equation to he integrated numerically, instead of 


(Y.3) • Tlie coefficients 


lP]f 


and 


can he calculated with much more 


accuracy than the factors involving 5, since they depend only on the 
two-hody solution. For analysis of error propagation, we write (Y.7) as 


?" = — [(p + 5) 

Pl 


£i 




-p] 


Ip 


M- 


F + ?' 

iP 


(Y.8) 
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The mechanics of the procedure, then, are easy to enumerate. The 
initial conditions are x© and .. q . Let 


p(to) = Xo 


xoUol 

p'(to) 

VM- 


Using these initial conditions, 
value of g to be considered. 


evaluate t. 



T7T 


p , p ' for each 


=* 56 * 0 . Using these initial conditions, integrate Equation 
(l*.?) to get §(g) and §'(g) . ?Tote that the first two terms on the 
right-hand side of Equation (Y.?) are functions of x and possibly x' . 
These are obtained by 

x(B) = p(e) + ?(b) 
x'(B) • p'(b) + :'(e) 

If at any point x is required, it can be found frcsn 


i[t(B)] = x'(B) 

|p(B)| 

Depending on the recltiflcatlon control logic, there will be places 
where the solution to Equation (Y.l) must be started over. At this 
point, the values t, x, x beconi, the new to, xq, and Xq, while g, f, 
and are reset to zero. 


(Y.Q) 


(y.io) 


(Y.ll) 



. u .'ison of Modes 

m r- — ■— ■ — I — — I 


a) It is inmediately apparent that eliminating the necessity of 
iteratively solving Equation (Y.2) will substantit Hj- increase tne 
speed of computation. 

b) An important advantage arises further from eliminating the some- 
times ponderous logic which supplies initial guesses for the 
iteratlv'i process and guarantees convergence of the solution. 

c) A third advantage of the 0-method is not quite so apparent, but no 
less in^>ortant. It is well known that the size of the integration 
time step can be increased as the distance from the center of 
attraction increases. Tbls change of the integration interval 
requires a cumbersome restart procedure . An examination of 
Equation (Y.5) shows that equal intervals of 0 correspond to 
tim ervals of increasing lengths as the distance incre--.ses . 

35ie time interval thus automatically expands and contracts correctly 
without outside intervention. 

d) Geometric stopping and printing conditions can ally be conveniently 
expressed in terms of 0, whereas they often require iterative deter- 
mination.' of the time. This advantage, however, is slight. 

e) If state vectors are required at fixed times, an iteration is neces- 
sary to find the corresponding value of 0 . In this case the 0- 
method is no better, ?=>nd ro worse, than the standard methods. If 
such vectors are required at frequent, closely spaced, time points 
(as in orbit determination, for instance), the advantage of the 
0-method is marginal. 
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N. SHADOW LOGIC 


A coordinate system is set up in the plane defined by the centers 
of the light-emitting source, the shadowing body, and the probe. Both 
bodies are assumed to be spherical, and hence all testing can be carried 
out in this plane. The diagram in Figure 3 shows this plane. 


Ihe coordinates are defined by unit vectors i and J: 

R . 

i = — ; 1*1 = 1; 1*1 = 0 (N.l) 

ct 


tj . - ^ ; J • i - 1 

where 


Vehicle coordinates in this system are given by: 


X 

V 


z 


V 


^c • K • 0 



(N.2) 


(N.3) 

(N.4) 

(H.5) 
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1 . Shadow Parameters 


a) The tl^iS of the umbra and penumbra cones are; 




cl 


^-1 

r 
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cl 
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b) The slopes of the bounding lines are: 


sin a = cos 0 » ~ 

u u d 
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■cos or = sin 9 
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f 12 


u=[-(^n 
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U cos (X 
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s a 9 


tan 0 = 
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u cos 9 
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c) Refraction Correction: (UMBRA.) 
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sin or = — T 
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sin 6p ■ cos cTp =j_l -{j-) J 


sin a 


tan or 


p cos (y 


sin 0 


teua 0 « 


p cos 0 


a' = a - € , 
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0' • 0 - e 
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sin or' = sin a cos e - cos a sin e 
u u u 


tan a ~ tan e 

tan or - ^ ^ tan e 

u u 


tan 0 - tan e 

®u ' 1 + tan 0 tan e 
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d) Befraction Correction: (PENUMBRA) 


Both € < 0 ?^ e > or_ j 
P P 


/y* a Of - e 
P P 


sin a' = I sin a cos e - cos a sin el 
pip P ' 


tan a ~ tan e 

tan or' » ^ — r 

T> 1 + tan a tan e 
P p 


tern 0 - tan e 

tan 0 * » 2 


P 1 + tan 0 tan e 
P 


d' = - slftn (tun or’) -—2 — r 
p p' sin a' 


The equations of the bounding lines are given below. 


2. The Testing Procedure 
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0 Sunlight 

3 Sunlight penumbra boundary 
0 Go to next test 


If 


Of exit here . 
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tan 9 ' 
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X 
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0 Penumbra 
< 0 Go to next test 
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> 0 Penumbra 

=s 0 Shadow penumbra boundary 
< 0 Shadow 


Qg and are stored and saved. The crossing times are fbund 
linearly interpolating for 0-values of Qg and respectively, to 
guarant -e that crossing from one region into another always occurs 
across these boundaries. 



0. SOLAR RADIATION PRESSURE 


The radiation pressure subroutine computes the force of solar 
radiation on the spacecraft if an appropriate pressure coefficient 
is used. The calculation relies on the shadow routine to set a 
trigger to multiply the pressure coefficient by 1.0, 0.5, or 0.0 for 
full sunlight, penumbra or umbra respectively. Therefore, the shadow 
subroutine must be used in conjunction with the radiation pressure 
routine for most cases. If the spacecraft is known to be continually 
in sunlight See Section VIII-2 to avoid elaborate shadow testing. 





m r' 


vs 


( 0 . 1 ) 


(See Section VIII-A-2 for definition of syaibols.) 


This radiation pressure subroutine has been found to be inexact 
for satellites of large area to mass ratio since it only controls the 
pressure to the nearest integration step. For such spacecraft (e.g. 
balloons) several degrees error in true ancanaly may result after 100 
days unless the integration is carried exactly to the boundaries. A 
modification to achieve this increased precision is available and will 
be included in future versions of the program. 
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p. 


ECLIPTIC COORDINATES 


Tbe ecliptic coordinates are an approximate «et obtained by a 
simple rotation of the eq.uatorial cooixiir' -^es about the x-axis through 
a fixed angle i = 23°26*35'' which io app: oximately the true obliquity 
for Jan 0.0, 1970. More exact coordinates may be obtained by changing 
SNE (unit normal to the ecliptic) as desired. 



Q. MOON ROTATING AND FIXED COORDINATE SYSTEM 


Geocentric coordinates of the vehicle based on the earth-moon plane 
are generated from the geocentric equatorial radius vector to the vehicle, 
R^, the geocentric unit vector in the direction of the moon and the 
vector in the direction of the moon's velocity, R^^. 


Coordinates in the rotating system, XROT etc., are found by using 
the current values of these vectors at each time step in the relations 


* >Ve • Bme 
^ ' ’Ve ■ (^ X 


(Q.l) 


where 



^ ^meI 


For the fixed axis system XINJ, etc., the initial vectors 
and Rjjg (to) ac the time of injection are used with the current 
value of Byj,. 
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R. TRAJECTORY SEARCH 


The prcjgram provides a search routine to obtain selected trajec- 
tories. The search is based on linear theory and varies the polar load 
input quanliities (independent variables) to search for desired dependent 
variables. A may^’"um of any six dependent variables may be selected. 
However, at this .Titing, the user is restricted to two dependent variable. 
Modifications to incorporate the remaining five are available and will 
be included in future versions of the program. The quantities at 

A A 

present are B • T and B • R. 


The components of the impact parameter vector (b • Tjjjp# B • Rj^jp) 
are referred to the ecliptic plane for T-planet trajectories and the 
Moon's orbital plane for lunar trajectories. The number of independent 
variables must equal the number of dependent variables for this routine 
to function. 

This version of the search routine is time consuming if the initial 
conditions are poorly approximated. Before using this routine, two 
things should be done. 



1. A first guess of the initial conditions of the nominal trajectory- 
should be obtained from a patched conic or a similar search 
program. 

2. The number of -variables should be kept to a minimum. It is lanned 
to automate the iteration scheme to go from two-body, to patched 
conic, to full trajectory, and to increase the number of variables 
to be adjusted, in optima.1 fashion. Even in its present form, 
however, it is extremely useful. 


The iterator uses a modified -version of the MIN-MAX Principle 
( Reference 6) . 


A. . 
ij 

Ax. 


^ii 


is the matrix of partials 

is the -vector of changes in the independent -variables 
is a diagonal matrix of weigh !:s 
is a -vector of residuals 


The system to be solved is 

^^Ji\j ^ii^^^i " '^ji^i 


^jAj ^ii ’ ®ij 


A 


ji^'i * h 
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Procedure 


The system 




is solved. If the value of is greater than SIZER, some arbitrary 
amount , set 





and solve the system again. Repeat these operations until is less 
than or equal to SIZER. Now, run a new nominal ctory with the new 

independent variable 


+ Ax^ 


a) If the new residuals y^ are greater than the previous ones, set 


^ ^ + X 


ij 


ij 


^ii 


Solve the system again and continue solving until the new residuals 
are less than the old. Now the system is re dy for a nev' iteration. 


b) If the new residuals y^ are less than the old, set 

Solve tie s^’tem again and continue solving until the new residue Is 
are greater or equal to the old. 


The iteration continues until eitner the maximum niunber of Iteratic ,iS 
(input) is exceeded or the residuals are less than or '^qual to an Inp’ t 
tolerance . 
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S. SaUATIOHS FOR FLIGHT PATH AZIMUTH AMD FLIGHT PATH AMGLE 


A subroutine computes the flight path azimuth and flight path suogle 
with the following equations; 

1. Flight path angle 


y ~ 



(S.l) 


^ is the vertical unit vector. In the geodetic system N is given by 


^ * [cos 9 cos (0 - 0o), cos 9 sin (9 - 0o, sir* 9 ] 


In the geocentric system cp is replaced by cp' . Alternatively, in the 
latter system 


A 

N 


R 

r 


2. Flight path azimuth 

A = sin"^ {v f (e-0o)}] 

(S.2) 

A 9 cos”^ I == i— - sin V sin cprl 

Lcos Y cos cp ^v ' 


Both formulas are used to determine the proper quadrant of A. To 
obtain the geocentric output, e^ • 0 , 9 is replaced by declination 
6 • 9' • 
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T. OSCUIA'JINfi ETJrMEWTS 


TLe osculating elements are obtained from the following equations: 


a = 

8 • f )" 

(T.l) 

n = 

CO 

1 

(T.2) 

e cos El /- r\ 

e cosh E J ’ V-'- ■ a^ 

(T.3) 

e sin El d 

e sinh Ei ~ f 

V |^ia| 

(T.4) 

M = 

r E - e sin E 
L e sinh E - E 

(T.5) 

t « 

P 

t 

n 

(T.6) 


The angles n, to, i are obtained from the vectors H auad P, where 


H » R X R 


(T.7) 


eP 



(T.8) 


In terms of these vectors. 


H 

cos i - in the first or fourth quadrant (T.9) 
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h sin 1 
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-H 
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h sin 1 


(T.IO) 


u) « P cos n + P sin n 
z y 

(T.n) 
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u. IMPACT PAKAMBTERS 


The "impact parameters" are coordinates in the "intact" plane. 

This pleme passes through the body (planet or the moon) and is normal to 
the Incoming asynqptote , The direction cosines of the asymptote are 
given by equations (U.l, U.2) in terms of unit vectors P (Appendix t) 
and 

(u.l) 
(U.2) 



In the plane defined by S as its normal, two unit vectors T.puD 

A A IMP 

^IMP defined. Tjjjp is parallel to the ecliptic plane for T-planet 
impacts, and to the mocn's orbital plane for moon iiQ)acts. Explicitly 


^IMP ^ 


A A 

H X S 

|S X s| 


(U.5) 


A 

where N is the unit normal to the ecliptic plane, or the moon's 

A A A 

orbital plane. is normal to both S and is the vector 

frcm the body to the vehicle as it crosses the impact plane. The data 
conputed are the dot products 


and 


®IMP ’ ^IMP 


^IMP * 
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V. MOON'S ORBITAL FLAME INPUT AMD OUTPUT 

A polar coordinate system is available for input and output which 
uses as its reference plEuie the moon's orbital plane and the vector 
from the moon to earth as unit '/ector. Polar coordinates in this 
system are defined analogous to geocentric polar coordinates. The 
cartesian coordinates in this system are computed by equations (H. 3 ) 
with 


and 


c = s = rp 
00 = 0 


Here r is the radius of the body of departure (earth or moon) . 

O 

These coordinates are then transformed to equatorial coordinates 
by a matrix C computed as follows; 

k ^ ;; 

3 = £ X £ 


The transformation matrix C is then given by 


C = 




(V.2) 


V-1 



and 




(V.J) 


Hie matrix C is unitary, and C ^ » C*, permitting easy inversion of 
equations (V.2) . 
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W. EQUATIONS FOR TRAUSLUNAR PLAHE INPUT 

The translunar plane input Is designed to permit easy visualization 
of the geometric relationships between initial conditions for circum- 
lunar trajectories and the motion of the moon. 

The initial conditions are given in a coordinate system referred to 
the translunar plane . This system has its x axis along the ascending 
node of the vehicle with respect to the moon's orbital plane, its y axis 
in the translunar plane at right angles to the ascending node, in the 
direction of motion. In this coordinate system, initial position and 
velocity vectors are given by 


XfL ’ h)cos Y 

7 tl =* (^B ^ 



= 0 


Here r is the radius of the body of departure (earth or moon) . 

B 

x^ » V sin (y - Y) 

» V cos (y - Y) (W.2) 



The translunar plane is positioned by giving its inclination i^ 
with respect to the moon's orbital plane and the lunar lead angle 0, 

the angle between the moon's position at injection and the descending 

» 

and may then be transformed into the equatorial 
system by the following series of rotations; 


node . The vectors 
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1. A rotation - 1^ about the axis will rotate the translunar plane 
into the moon's orbital plane. 

2. A rotation of n - + cp) about the new z-axis will refer the moon's 

orbital plane coordinate systen to the ascending node of the moon's orbital 
plane (with respect to the equator) as x-axis . 

Here X^, stands for the argument of latitude of the moon. These 
M 

rotations are performed by multiplying and 


by the matrix: 


■cos + cp) 
•sin (Xj^ + cp) 
0 


sin + cp) 

•cos cp)cos i^ 


■sin «p) sin 1^^^ 

cos (Xj^ ■^ cp)sin i^ 


sin 1 


TL 


cos i 


TL 


(W.5) 


5- The moon's orbital plane (MOP) is rotated about its node through an 
angle - 1„ (the inclination of the MOP) 

!<-. The ascending node is brought into coincidence with the vernal 
equinox by a rotation - These two rotations are embodied in the 
matrix 


B 


cos 

sinP^ 

0 


and thus: 


cos 
sin 1 


cos 

cos 
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sin fij, sin i^ 
-cos sin i^ 
cos 1 
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R . (BA) 
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